]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Added the possibility to compute laplacians of a FEFieldFunction.
authorgoll <goll@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Feb 2012 15:16:54 +0000 (15:16 +0000)
committergoll <goll@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 16 Feb 2012 15:16:54 +0000 (15:16 +0000)
git-svn-id: https://svn.dealii.org/trunk@25100 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/include/deal.II/numerics/fe_field_function.h
deal.II/include/deal.II/numerics/fe_field_function.templates.h

index 151c5255e9163df3f8bc3a6547b5fe0a789e56bb..fcb22bdfee30b4da8052cba4465a8f9887651f2c 100644 (file)
@@ -297,6 +297,41 @@ namespace Functions
       gradient_list (const std::vector< Point< dim > > &p, 
                     std::vector< Tensor< 1, dim > > &gradients, 
                     const unsigned int component=0) const;
+
+
+      /**
+       * Compute the Laplacian of a
+       * given component at point <tt>p</tt>.
+       */
+      virtual double
+      laplacian (const Point<dim>   &p,
+          const unsigned int  component = 0) const;
+
+      /**
+       * Compute the Laplacian of all
+       * components at point <tt>p</tt> and
+       * store them in <tt>values</tt>.
+       */
+      virtual void
+      vector_laplacian (const Point<dim>   &p,
+          Vector<double>     &values) const;
+
+      /**
+       * Compute the Laplacian of one
+       * component at a set of points.
+       */
+      virtual void
+      laplacian_list (const std::vector<Point<dim> > &points,
+          std::vector<double>            &values,
+          const unsigned int              component = 0) const;
+
+      /**
+       * Compute the Laplacians of all
+       * components at a set of points.
+       */
+      virtual void
+      vector_laplacian_list (const std::vector<Point<dim> > &points,
+          std::vector<Vector<double> >   &values) const;
   
                                       /**
                                        * Create quadrature
index 70dce9ed4417cf7efc9c41d5f56cbcf2303de667..a5687b4e289e831abe9f2f0496b058a12f2c34a8 100644 (file)
@@ -139,6 +139,52 @@ namespace Functions
     return grads[comp];
   }
 
+
+
+  template <int dim, typename DH, typename VECTOR>
+  void FEFieldFunction<dim, DH, VECTOR>::vector_laplacian
+  (const Point<dim> &p,Vector<double>   &values) const
+  {
+    Assert (values.size() == n_components,
+        ExcDimensionMismatch(values.size(), n_components));
+    typename DH::active_cell_iterator cell = cell_hint.get();
+    if (cell == dh->end())
+      cell = dh->begin_active();
+    Point<dim> qp = mapping.transform_real_to_unit_cell(cell, p);
+
+    // Check if we already have all we need
+    if (!GeometryInfo<dim>::is_inside_unit_cell(qp))
+      {
+        const std::pair<typename DH::active_cell_iterator, Point<dim> > my_pair
+        = GridTools::find_active_cell_around_point (mapping, *dh, p);
+        cell = my_pair.first;
+        qp = my_pair.second;
+      }
+
+    cell_hint.get() = cell;
+
+    // Now we can find out about the point
+    Quadrature<dim> quad(qp);
+    FEValues<dim> fe_v(mapping, dh->get_fe(), quad,
+        update_hessians);
+    fe_v.reinit(cell);
+    std::vector< Vector<double> > vvalues (1, values);
+    fe_v.get_function_laplacians(data_vector, vvalues);
+    values = vvalues[0];
+  }
+
+
+
+  template <int dim, typename DH, typename VECTOR>
+  double FEFieldFunction<dim, DH, VECTOR>::laplacian
+  (const Point<dim>   &p, const unsigned int  comp) const
+  {
+    Vector<double> lap(n_components);
+    vector_laplacian(p, lap);
+    return lap[comp];
+  }
+
+
                                   // Now the list versions
                                   // ==============================
 
@@ -249,6 +295,58 @@ namespace Functions
       values[q] = vvalues[q][component];
   }
 
+
+  template <int dim, typename DH, typename VECTOR>
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  vector_laplacian_list (const std::vector<Point< dim > > &    points,
+         std::vector< Vector<double> > &values) const
+  {
+    Assert(points.size() == values.size(),
+     ExcDimensionMismatch(points.size(), values.size()));
+
+    std::vector<typename DH::active_cell_iterator > cells;
+    std::vector<std::vector<Point<dim> > > qpoints;
+    std::vector<std::vector<unsigned int> > maps;
+
+    unsigned int ncells = compute_point_locations(points, cells, qpoints, maps);
+
+             // Now gather all the informations we need
+    for (unsigned int i=0; i<ncells; ++i)
+      {
+           // Number of quadrature points on this cell
+  unsigned int nq = qpoints[i].size();
+
+           // Construct a quadrature formula
+  std::vector< double > ww(nq, 1./((double) nq));
+  Quadrature<dim> quad(qpoints[i], ww);
+
+           // Get a function value object
+  FEValues<dim> fe_v(mapping, dh->get_fe(), quad,
+         update_hessians);
+  fe_v.reinit(cells[i]);
+  std::vector< Vector<double> > vvalues (nq, Vector<double>(n_components));
+  fe_v.get_function_laplacians(data_vector, vvalues);
+  for (unsigned int q=0; q<nq; ++q)
+    values[maps[i][q]] = vvalues[q];
+      }
+  }
+
+  template <int dim, typename DH, typename VECTOR>
+  void
+  FEFieldFunction<dim, DH, VECTOR>::
+  laplacian_list (const std::vector<Point< dim > > &points,
+        std::vector< double > &values,
+        const unsigned int  component) const
+  {
+    Assert(points.size() == values.size(),
+     ExcDimensionMismatch(points.size(), values.size()));
+    std::vector< Vector<double> > vvalues(points.size(), Vector<double>(n_components));
+    vector_laplacian_list(points, vvalues);
+    for (unsigned int q=0; q<points.size(); ++q)
+      values[q] = vvalues[q](component);
+  }
+
   
 
   template <int dim, typename DH, typename VECTOR>

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.