]> https://gitweb.dealii.org/ - dealii.git/commitdiff
More numerics and more efficient version of FEValues::reinit
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 1 Apr 1998 12:37:54 +0000 (12:37 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 1 Apr 1998 12:37:54 +0000 (12:37 +0000)
git-svn-id: https://svn.dealii.org/trunk@108 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/poisson/poisson.cc
deal.II/deal.II/include/fe/fe.h
deal.II/deal.II/include/fe/fe_lib.lagrange.h
deal.II/deal.II/include/numerics/assembler.h
deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/source/fe/fe.cc
deal.II/deal.II/source/fe/fe_lib.linear.cc
deal.II/deal.II/source/numerics/assembler.cc
deal.II/deal.II/source/numerics/base.cc
tests/big-tests/poisson/poisson.cc

index a1f644440822e5915d8d0670e5f213ecb7061aef..03fa5a31ab8f200c63c96349235fb664aa03f713 100644 (file)
@@ -131,14 +131,14 @@ int main () {
 //  tria.create_hyper_ball(Point<2>(2,3),4);
 //  tria.set_boundary (&boundary);
   
-/*  tria.refine_global (1);
-  (--tria.last_active())->set_refine_flag();
-  tria.execute_refinement ();
-  tria.begin_active(2)->set_refine_flag();
-  tria.execute_refinement ();
+  tria.refine_global (1);
+//  (--tria.last_active())->set_refine_flag();
+//  tria.execute_refinement ();
+//  tria.begin_active(2)->set_refine_flag();
+//  tria.execute_refinement ();
   tria.refine_global (2);
-*/
 
+/*
   const unsigned int dim=2;
   tria.refine_global (1);
        
@@ -163,7 +163,7 @@ int main () {
       tria.execute_refinement ();
     };
   tria.refine_global (1);
-  
+*/
 
   cout << "Distributing dofs... "; 
   dof.distribute_dofs (fe);
index f74d3e841400d8de9b404f6b142266f9dcfffefc..f3ffc7d85cdbb47daf014c4b8a61b74c501bbbe5 100644 (file)
@@ -41,10 +41,70 @@ template <int dim> class Quadrature;
   $dx dy$ from the real to the unit cell, we have to take the determinant of
   the inverse matrix, which is the reciprocal value of the determinant of the
   matrix defined above.
+
+  The #FEValues# object keeps track of those fields which really need to
+  be computed, since the computation of the gradients of the ansatz functions
+  on each real cell can be quite an expensive thing if it is not needed. The
+  object knows about which fields are needed by the #UpdateStruct# object
+  passed through the constructor. In debug mode, the accessor functions, which
+  return values from the different fields, check whether the required field
+  was initialized, thus avoiding use of unitialized data.
   */
 template <int dim>
 class FEValues {
   public:
+                                    /**
+                                     * Provide a structure which tells the
+                                     * #reinit# function, which fields are
+                                     * to be updated for each cell. E.g. if
+                                     * you do not need the gradients since
+                                     * you want to assemble the mass matrix,
+                                     * you can switch that off. By default,
+                                     * all flags are off, i.e. no
+                                     * reinitialization will be done.
+                                     *
+                                     * A structure of this type has to be
+                                     * passed to the constructor of the
+                                     * #FEValues# object. 
+                                     */
+    struct UpdateStruct {
+                                        /**
+                                         * Constructor. Sets all fields to
+                                         * false.
+                                         */
+       UpdateStruct ();
+                                        /**
+                                         * Compute quadrature points in real
+                                         * space (not on unit cell).
+                                         */
+       bool update_q_points;
+                                        /**
+                                         * Transform gradients on unit cell to
+                                         * gradients on real cell.
+                                         */
+       bool update_gradients;
+                                        /**
+                                         * Compute jacobian matrices of the
+                                         * transform between unit and real cell
+                                         * in the evaluation points.
+                                         */
+       bool update_jacobians;
+                                        /**
+                                         * Compute the JxW values (Jacobian
+                                         * determinant at the quadrature point
+                                         * times the weight of this point).
+                                         */
+       bool update_JxW_values;
+                                        /**
+                                         * Compute the points on the real cell
+                                         * on which the ansatz functions are
+                                         * located.
+                                         */
+       bool update_ansatz_points;
+    };
+
+    
+    
                                     /**
                                      * Number of quadrature points.
                                      */
@@ -63,7 +123,8 @@ class FEValues {
                                      * quadrature rule.
                                      */
     FEValues (const FiniteElement<dim> &,
-             const Quadrature<dim> &);
+             const Quadrature<dim> &,
+             const UpdateStruct &);
 
                                     /**
                                      * Return the value of the #i#th shape
@@ -147,6 +208,14 @@ class FEValues {
                    int, int,
                    << "The index " << arg1
                    << " is out of range, it should be less than " << arg2);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcAccessToUninitializedField);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcCannotInitializeField);
     
   private:
                                     /**
@@ -226,6 +295,12 @@ class FEValues {
                                      * is set each time #reinit# is called.
                                      */
     vector<dFMatrix>     jacobi_matrices;
+
+                                    /**
+                                     * Store which fields are to be updated by
+                                     * the reinit function.
+                                     */
+    UpdateStruct         update_flags;
 };
 
 
@@ -344,6 +419,14 @@ class FiniteElementBase {
                                      * elements need different transformations
                                      * of the unit cell to a real cell.
                                      *
+                                     * The computation of the two fields may
+                                     * share some common code, which is why we
+                                     * put it in one function. However, it may
+                                     * not always be necessary to really
+                                     * compute both fields, so there are two
+                                     * bool flags which tell the function which
+                                     * of the fields to actually compute.
+                                     *
                                      * Refer to the documentation of the
                                      * \Ref{FEValues} class for a definition
                                      * of the Jacobi matrix.
@@ -357,7 +440,9 @@ class FiniteElementBase {
     virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
-                                vector<Point<dim> > &points) const;
+                                const bool           compute_jacobians,
+                                vector<Point<dim> > &points,
+                                const bool           compute_points) const;
     
                                     /**
                                      * Comparison operator. We also check for
@@ -547,7 +632,9 @@ class FiniteElement<1> : public FiniteElementBase<1> {
     virtual void fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                 const vector<Point<1> >               &unit_points,
                                 vector<dFMatrix>  &jacobians,
-                                vector<Point<1> > &points) const;
+                                const bool         compute_jacobians,
+                                vector<Point<1> > &points,
+                                const bool         compute_points) const;
 };
 
 
@@ -667,7 +754,9 @@ class FiniteElement<2> : public FiniteElementBase<2> {
     virtual void fill_fe_values (const Triangulation<2>::cell_iterator &cell,
                                 const vector<Point<2> >               &unit_points,
                                 vector<dFMatrix>  &jacobians,
-                                vector<Point<2> > &points) const;
+                                const bool         compute_jacobians,
+                                vector<Point<2> > &points,
+                                const bool         compute_points) const;
 };
 
 
@@ -691,6 +780,7 @@ template <int dim>
 inline
 const vector<vector<Point<dim> > > &
 FEValues<dim>::get_shape_grads () const {
+  Assert (update_flags.update_gradients, ExcAccessToUninitializedField());
   return shape_gradients;
 };
 
@@ -700,6 +790,7 @@ template <int dim>
 inline
 const vector<Point<dim> > &
 FEValues<dim>::get_quadrature_points () const {
+  Assert (update_flags.update_q_points, ExcAccessToUninitializedField());
   return quadrature_points;
 };
 
@@ -709,6 +800,7 @@ template <int dim>
 inline
 const vector<double> &
 FEValues<dim>::get_JxW_values () const {
+  Assert (update_flags.update_JxW_values, ExcAccessToUninitializedField());
   return JxW_values;
 };
 
index e14f219e403f194fac65527726193356fc03daab..279f872740ad03f8e343380bd81443ba7e4bb08b 100644 (file)
@@ -70,7 +70,9 @@ class FELinear : public FiniteElement<dim> {
     virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
-                                vector<Point<dim> > &points) const;
+                                const bool           compute_jacobians,
+                                vector<Point<dim> > &points,
+                                const bool           compute_points) const;
 };
 
 
@@ -129,7 +131,9 @@ class FEQuadratic : public FiniteElement<dim> {
     virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
-                                vector<Point<dim> > &points) const;
+                                const bool           compute_jacobians,
+                                vector<Point<dim> > &points,
+                                const bool           compute_points) const;
 };
 
 
@@ -188,7 +192,9 @@ class FECubic : public FiniteElement<dim> {
     virtual void fill_fe_values (const Triangulation<dim>::cell_iterator &cell,
                                 const vector<Point<dim> >               &unit_points,
                                 vector<dFMatrix>    &jacobians,
-                                vector<Point<dim> > &points) const;
+                                const bool           compute_jacobians,
+                                vector<Point<dim> > &points,
+                                const bool           compute_points) const;
 };
 
 
index bd2192629f383d4b5df6f274de8394fcdcebfbc3..44af7879eb023bec15cccfac147573852ba2551a 100644 (file)
@@ -105,7 +105,8 @@ struct AssemblerData {
                   const bool                assemble_rhs,
                   ProblemBase<dim>         &problem,
                   const Quadrature<dim>    &quadrature,
-                  const FiniteElement<dim> &fe);
+                  const FiniteElement<dim> &fe,
+                  const FEValues<dim>::UpdateStruct &update_flags);
     
                                     /**
                                      * Pointer to the dof handler object
@@ -135,6 +136,12 @@ struct AssemblerData {
                                      * object.
                                      */
     const FiniteElement<dim> &fe;
+                                    /**
+                                     * Store which of the fields of the
+                                     * FEValues object need to be reinitialized
+                                     * on each cell.
+                                     */
+    const FEValues<dim>::UpdateStruct update_flags;
 };
 
 
index 572239610dcba04f289b3c213509cb341c3408c5..cfc2d4297db730f5ea28b29864599bc6d416fe3a 100644 (file)
@@ -248,10 +248,11 @@ class ProblemBase {
                                      * For what exactly happens here, refer to
                                      * the general doc of this class.
                                      */
-    virtual void assemble (const Equation<dim>      &equation,
-                          const Quadrature<dim>    &q,
-                          const FiniteElement<dim> &fe,
-                          const DirichletBC        &dirichlet_bc = DirichletBC());
+    virtual void assemble (const Equation<dim>               &equation,
+                          const Quadrature<dim>             &q,
+                          const FiniteElement<dim>          &fe,
+                          const FEValues<dim>::UpdateStruct &update_flags,
+                          const DirichletBC                 &dirichlet_bc = DirichletBC());
     
                                     /**
                                      * Solve the system of equations.
index 7c591c7c4e4528c528004146d040ec1a4070798c..8e0895ebc073db3179313269d7db00371545bee3 100644 (file)
@@ -7,12 +7,25 @@
 
 
 
+template <int dim>
+FEValues<dim>::UpdateStruct::UpdateStruct () :
+               update_q_points(false),
+               update_gradients(false),
+               update_jacobians(false),
+               update_JxW_values(false),
+               update_ansatz_points(false) {};
+               
+
+
+
+
 /*------------------------------- FEValues -------------------------------*/
 
 
 template <int dim>
 FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
-                        const Quadrature<dim>    &quadrature) :
+                        const Quadrature<dim>    &quadrature,
+                        const UpdateStruct       &update_flags) :
                n_quadrature_points(quadrature.n_quadrature_points),
                total_dofs(fe.total_dofs),
                shape_values(fe.total_dofs, quadrature.n_quadrature_points),
@@ -25,7 +38,8 @@ FEValues<dim>::FEValues (const FiniteElement<dim> &fe,
                quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
                unit_quadrature_points(quadrature.n_quadrature_points, Point<dim>()),
                jacobi_matrices (quadrature.n_quadrature_points,
-                                dFMatrix(dim,dim))
+                                dFMatrix(dim,dim)),
+               update_flags (update_flags)
 {
   for (unsigned int i=0; i<fe.total_dofs; ++i)
     for (unsigned int j=0; j<n_quadrature_points; ++j) 
@@ -61,6 +75,7 @@ FEValues<dim>::shape_grad (const unsigned int i,
                           const unsigned int j) const {
   Assert (i<(unsigned int)shape_values.m(), ExcInvalidIndex (i, shape_values.m()));
   Assert (j<(unsigned int)shape_values.n(), ExcInvalidIndex (j, shape_values.n()));
+  Assert (update_flags.update_gradients, ExcAccessToUninitializedField());
 
   return shape_gradients[i][j];
 };
@@ -70,6 +85,7 @@ FEValues<dim>::shape_grad (const unsigned int i,
 template <int dim>
 const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+  Assert (update_flags.update_q_points, ExcAccessToUninitializedField());
   
   return quadrature_points[i];
 };
@@ -79,6 +95,7 @@ const Point<dim> & FEValues<dim>::quadrature_point (const unsigned int i) const
 template <int dim>
 double FEValues<dim>::JxW (const unsigned int i) const {
   Assert (i<n_quadrature_points, ExcInvalidIndex(i, n_quadrature_points));
+  Assert (update_flags.update_JxW_values, ExcAccessToUninitializedField());
   
   return JxW_values[i];
 };
@@ -89,35 +106,50 @@ template <int dim>
 void FEValues<dim>::reinit (const typename Triangulation<dim>::cell_iterator &cell,
                            const FiniteElement<dim>                         &fe) {
                                   // fill jacobi matrices and real
-                                  //quadrature points
-  fe.fill_fe_values (cell,
-                    unit_quadrature_points,
-                    jacobi_matrices,
-                    quadrature_points);
-
-                                  // compute gradients on real element
-  for (unsigned int i=0; i<fe.total_dofs; ++i)
-    for (unsigned int j=0; j<n_quadrature_points; ++j)
-      for (unsigned int s=0; s<dim; ++s)
-       {
-         shape_gradients[i][j](s) = 0;
-
-                                          // (grad psi)_s =
-                                          // (grad_{\xi\eta})_b J_{bs}
-                                          // with J_{bs}=(d\xi_b)/(dx_s)
-         for (unsigned int b=0; b<dim; ++b)
-           shape_gradients[i][j](s)
-             +=
-             unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
-       };
+                                  // quadrature points
+  if (update_flags.update_jacobians || update_flags.update_q_points)
+    fe.fill_fe_values (cell,
+                      unit_quadrature_points,
+                      jacobi_matrices,
+                      update_flags.update_jacobians,
+                      quadrature_points,
+                      update_flags.update_q_points);
+
+                                  // compute gradients on real element if
+                                  // requested
+  if (update_flags.update_gradients) 
+    {
+      Assert (update_flags.update_jacobians, ExcCannotInitializeField());
+      
+      for (unsigned int i=0; i<fe.total_dofs; ++i)
+       for (unsigned int j=0; j<n_quadrature_points; ++j)
+         for (unsigned int s=0; s<dim; ++s)
+           {
+             shape_gradients[i][j](s) = 0;
+             
+                                              // (grad psi)_s =
+                                              // (grad_{\xi\eta})_b J_{bs}
+                                              // with J_{bs}=(d\xi_b)/(dx_s)
+             for (unsigned int b=0; b<dim; ++b)
+               shape_gradients[i][j](s)
+                 +=
+                 unit_shape_gradients[i][j](b) * jacobi_matrices[j](b,s);
+           };
+    };
+  
   
                                   // compute Jacobi determinants in
                                   // quadrature points.
                                   // refer to the general doc for
                                   // why we take the inverse of the
                                   // determinant
-  for (unsigned int i=0; i<n_quadrature_points; ++i)
-    JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
+  if (update_flags.update_JxW_values) 
+    {
+      Assert (update_flags.update_jacobians,
+             ExcCannotInitializeField());
+      for (unsigned int i=0; i<n_quadrature_points; ++i)
+       JxW_values[i] = weights[i] / jacobi_matrices[i].determinant();
+    };
 };
 
 
@@ -199,7 +231,9 @@ template <int dim>
 void FiniteElementBase<dim>::fill_fe_values (const typename Triangulation<dim>::cell_iterator &,
                                             const vector<Point<dim> > &,
                                             vector<dFMatrix> &,
-                                            vector<Point<dim> > &) const {
+                                            const bool,
+                                            vector<Point<dim> > &,
+                                            const bool) const {
   Assert (false, ExcPureFunctionCalled());
 };
 
@@ -219,14 +253,18 @@ bool FiniteElement<1>::operator == (const FiniteElement<1> &f) const {
 void FiniteElement<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                       const vector<Point<1> > &unit_points,
                                       vector<dFMatrix>  &jacobians,
-                                      vector<Point<1> > &points) const {
+                                      const bool         compute_jacobians,
+                                      vector<Point<1> > &points,
+                                      const bool         compute_points) const {
                                   // local mesh width
   const double h=(cell->vertex(1)(0) - cell->vertex(0)(0));
 
   for (unsigned int i=0; i<points.size(); ++i) 
     {
-      jacobians[i](0,0) = 1./h;
-      points[i] = cell->vertex(0) + h*unit_points[i];
+      if (compute_jacobians)
+       jacobians[i](0,0) = 1./h;
+      if (compute_points)
+       points[i] = cell->vertex(0) + h*unit_points[i];
     };
 };
 
@@ -244,7 +282,9 @@ bool FiniteElement<2>::operator == (const FiniteElement<2> &f) const {
 void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
                                       const vector<Point<2> > &,
                                       vector<dFMatrix> &,
-                                      vector<Point<2> > &) const {
+                                      const bool,
+                                      vector<Point<2> > &,
+                                      const bool) const {
   Assert (false, ExcPureFunctionCalled());
 };
 
@@ -254,6 +294,9 @@ void FiniteElement<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
 
 /*------------------------------- Explicit Instantiations -------------*/
 
+template struct FEValues<1>::UpdateStruct;
+template struct FEValues<2>::UpdateStruct;
+
 template class FEValues<1>;
 template class FEValues<2>;
 
index d0a167698b350a1bbc059dfb3670ba6122b69809..a9c6fd16a27668f8f97e652336b490e14e9fc008 100644 (file)
@@ -79,9 +79,13 @@ FELinear<1>::shape_grad(const unsigned int i,
 void FELinear<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                  const vector<Point<1> >               &unit_points,
                                  vector<dFMatrix>  &jacobians,
-                                 vector<Point<1> > &points) const {
+                                 const bool         compute_jacobians,
+                                 vector<Point<1> > &points,
+                                 const bool         compute_points) const {
                                   // simply pass down
-  FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+  FiniteElement<1>::fill_fe_values (cell, unit_points,
+                                   jacobians, compute_jacobians,
+                                   points, compute_points);
 };
 
 
@@ -227,32 +231,37 @@ FELinear<2>::shape_grad (const unsigned int i,
 void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
                                  const vector<Point<2> >               &unit_points,
                                  vector<dFMatrix>  &jacobians,
-                                 vector<Point<2> > &points) const {
+                                 const bool         compute_jacobians,
+                                 vector<Point<2> > &points,
+                                 const bool         compute_points) const {
   const unsigned int dim=2;
   const unsigned int n_vertices=4;
   unsigned int n_points=unit_points.size();
 
-                                  // recomputation of values
   Point<dim> vertices[n_vertices];
   for (unsigned int l=0; l<n_vertices; ++l)
     vertices[l] = cell->vertex(l);
   
+
+  if (compute_points) 
+    {
+                                      // initialize points to zero
+      for (unsigned int i=0; i<n_points; ++i)
+       points[i] = Point<dim> ();
+      
+                                      // note: let x_l be the vector of the
+                                      // lth quadrature point in real space and
+                                      // xi_l that on the unit cell, let further
+                                      // p_j be the vector of the jth vertex
+                                      // of the cell in real space and
+                                      // N_j(xi_l) be the value of the associated
+                                      // basis function at xi_l, then
+                                      // x_l(xi_l) = sum_j p_j N_j(xi_l)
+      for (unsigned int j=0; j<n_vertices; ++j) 
+       for (unsigned int l=0; l<n_points; ++l) 
+         points[l] += vertices[j] * shape_value(j, unit_points[l]);
+    };
   
-                                  // initialize points to zero
-  for (unsigned int i=0; i<n_points; ++i)
-    points[i] = Point<dim> ();
-  
-                                  // note: let x_l be the vector of the
-                                  // lth quadrature point in real space and
-                                  // xi_l that on the unit cell, let further
-                                  // p_j be the vector of the jth vertex
-                                  // of the cell in real space and
-                                  // N_j(xi_l) be the value of the associated
-                                  // basis function at xi_l, then
-                                  // x_l(xi_l) = sum_j p_j N_j(xi_l)
-  for (unsigned int j=0; j<n_vertices; ++j) 
-    for (unsigned int l=0; l<n_points; ++l) 
-      points[l] += vertices[j] * shape_value(j, unit_points[l]);
 
 /* jacobi matrices: compute d(x)/d(xi) and invert this
    Let M(l) be the inverse of J at the quadrature point l, then
@@ -272,24 +281,27 @@ void FELinear<2>::fill_fe_values (const Triangulation<2>::cell_iterator &cell,
   However, we rewrite the loops to only compute the gradient once for
   each integration point and basis function.
 */
-  dFMatrix M(dim,dim);
-  for (unsigned int l=0; l<n_points; ++l) 
+  if (compute_jacobians) 
     {
-      M.clear ();
-      for (unsigned int s=0; s<n_vertices; ++s)
+      dFMatrix M(dim,dim);
+      for (unsigned int l=0; l<n_points; ++l) 
        {
-                                          // we want a linear transform and
-                                          // if we prepend the class name in
-                                          // front of the #shape_grad#, we
-                                          // need not use virtual function
-                                          // calls.
-         const Point<dim> gradient
-           = FELinear<dim>::shape_grad (s, unit_points[l]);
-         for (unsigned int i=0; i<dim; ++i)
-           for (unsigned int j=0; j<dim; ++j)
-             M(i,j) += vertices[s](i) * gradient(j);
+         M.clear ();
+         for (unsigned int s=0; s<n_vertices; ++s)
+           {
+                                              // we want a linear transform and
+                                              // if we prepend the class name in
+                                              // front of the #shape_grad#, we
+                                              // need not use virtual function
+                                              // calls.
+             const Point<dim> gradient
+               = FELinear<dim>::shape_grad (s, unit_points[l]);
+             for (unsigned int i=0; i<dim; ++i)
+               for (unsigned int j=0; j<dim; ++j)
+                 M(i,j) += vertices[s](i) * gradient(j);
+           };
+         jacobians[l].invert(M);
        };
-      jacobians[l].invert(M);
     };
 };
 
@@ -306,9 +318,13 @@ FEQuadratic<1>::FEQuadratic () :
 void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                     const vector<Point<1> >               &unit_points,
                                     vector<dFMatrix>  &jacobians,
-                                    vector<Point<1> > &points) const {
+                                    const bool         compute_jacobians,
+                                    vector<Point<1> > &points,
+                                    const bool         compute_points) const {
                                   // simply pass down
-  FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+  FiniteElement<1>::fill_fe_values (cell, unit_points,
+                                   jacobians, compute_jacobians,
+                                   points, compute_points);
 };
 
 
@@ -355,7 +371,9 @@ FEQuadratic<dim>::shape_grad (const unsigned int i,
 void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
                                     const vector<Point<2> >               &,
                                     vector<dFMatrix>  &,
-                                    vector<Point<2> > &) const {
+                                    const bool,
+                                    vector<Point<2> > &,
+                                    const bool) const {
   Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
 };
 
@@ -372,9 +390,13 @@ FECubic<1>::FECubic () :
 void FECubic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell,
                                 const vector<Point<1> >               &unit_points,
                                 vector<dFMatrix>  &jacobians,
-                                vector<Point<1> > &points) const {
+                                const bool         compute_jacobians,
+                                vector<Point<1> > &points,
+                                const bool          compute_points) const {
                                   // simply pass down
-  FiniteElement<1>::fill_fe_values (cell, unit_points, jacobians, points);
+  FiniteElement<1>::fill_fe_values (cell, unit_points,
+                                   jacobians, compute_jacobians,
+                                   points, compute_points);
 };
 
 
@@ -413,7 +435,9 @@ FECubic<dim>::shape_grad (const unsigned int i,
 void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &,
                                 const vector<Point<2> >               &,
                                 vector<dFMatrix>  &,
-                                vector<Point<2> > &) const {
+                                const bool,
+                                vector<Point<2> > &,
+                                const bool) const {
   Assert (false, typename FiniteElementBase<2>::ExcNotImplemented());
 };
 
index 218a66cddb5a8e5786d0c05d5c6598baf3d30ce8..591e2de10d1f3c8ab519d4ef53d439e1180bc0d8 100644 (file)
@@ -50,13 +50,15 @@ AssemblerData<dim>::AssemblerData (const DoFHandler<dim>   &dof,
                                   const bool               assemble_rhs,
                                   ProblemBase<dim>         &problem,
                                   const Quadrature<dim>    &quadrature,
-                                  const FiniteElement<dim> &fe) :
+                                  const FiniteElement<dim> &fe,
+                                  const FEValues<dim>::UpdateStruct &update_flags) :
                dof(dof),
                assemble_matrix(assemble_matrix),
                assemble_rhs(assemble_rhs),
                problem(problem),
                quadrature(quadrature),
-               fe(fe) {};
+               fe(fe),
+               update_flags(update_flags) {};
 
 
 
@@ -75,7 +77,8 @@ Assembler<dim>::Assembler (Triangulation<dim> *tria,
                problem (((AssemblerData<dim>*)local_data)->problem),
                fe(((AssemblerData<dim>*)local_data)->fe),
                fe_values (((AssemblerData<dim>*)local_data)->fe,
-                          ((AssemblerData<dim>*)local_data)->quadrature)
+                          ((AssemblerData<dim>*)local_data)->quadrature,
+                          ((AssemblerData<dim>*)local_data)->update_flags)
 {
   Assert ((unsigned int)problem.system_matrix.m() == dof_handler->n_dofs(),
          ExcInvalidData());
index 84c946b70867d3d13f0db234da3dab48fb71461f..64bc08bfc7404e4ad6a8eeba218cd26007e9a3d2 100644 (file)
@@ -45,10 +45,11 @@ ProblemBase<dim>::~ProblemBase () {};
 
 
 template <int dim>
-void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
-                                const Quadrature<dim>    &quadrature,
-                                const FiniteElement<dim> &fe,
-                                const DirichletBC        &dirichlet_bc) {
+void ProblemBase<dim>::assemble (const Equation<dim>               &equation,
+                                const Quadrature<dim>             &quadrature,
+                                const FiniteElement<dim>          &fe,
+                                const FEValues<dim>::UpdateStruct &update_flags,
+                                const DirichletBC                 &dirichlet_bc) {
   system_sparsity.reinit (dof_handler->n_dofs(),
                          dof_handler->n_dofs(),
                          dof_handler->max_couplings_between_dofs());
@@ -66,18 +67,18 @@ void ProblemBase<dim>::assemble (const Equation<dim>      &equation,
                                   // reinit solution vector, preset
                                   // with zeroes.
   solution.reinit (dof_handler->n_dofs());
-
+  
                                   // create assembler
   AssemblerData<dim> data (*dof_handler,
                           true, true, //assemble matrix and rhs
                           *this,
                           quadrature,
-                          fe);
+                          fe,
+                          update_flags);
   TriaActiveIterator<dim, Assembler<dim> > assembler (tria,
                                                      tria->begin_active()->level(),
                                                      tria->begin_active()->index(),
                                                      &data);
-
                                   // loop over all cells, fill matrix and rhs
   do 
     {
@@ -128,7 +129,10 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
   difference.erase (difference.begin(), difference.end());
   difference.reserve (tria->n_cells());
 
-  FEValues<dim> fe_values(fe, q);
+  FEValues<dim>::UpdateStruct update_flags;
+  update_flags.update_q_points   = true;
+  update_flags.update_JxW_values = true;
+  FEValues<dim> fe_values(fe, q, update_flags);
   
                                   // loop over all cells
                                   // (we need an iterator on the triangulation for
@@ -187,9 +191,10 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &exact_sol
                                       psi);
                                             // then subtract finite element
                                             // solution
-           for (unsigned int j=0; j<n_dofs; ++j)
+           const vector<double> &JxW_values = fe_values.get_JxW_values();
+           for (unsigned int j=0; j<n_dofs; ++j) 
              for (unsigned int i=0; i<n_dofs; ++i)
-               psi[j] -= dof_values[i]*shape_values(i,j);
+               psi[j] -= dof_values[i]*shape_values(i,j)*JxW_values[j];
 
                                             // for L1_norm and Linfty_norm:
                                             // take absolute
index a1f644440822e5915d8d0670e5f213ecb7061aef..03fa5a31ab8f200c63c96349235fb664aa03f713 100644 (file)
@@ -131,14 +131,14 @@ int main () {
 //  tria.create_hyper_ball(Point<2>(2,3),4);
 //  tria.set_boundary (&boundary);
   
-/*  tria.refine_global (1);
-  (--tria.last_active())->set_refine_flag();
-  tria.execute_refinement ();
-  tria.begin_active(2)->set_refine_flag();
-  tria.execute_refinement ();
+  tria.refine_global (1);
+//  (--tria.last_active())->set_refine_flag();
+//  tria.execute_refinement ();
+//  tria.begin_active(2)->set_refine_flag();
+//  tria.execute_refinement ();
   tria.refine_global (2);
-*/
 
+/*
   const unsigned int dim=2;
   tria.refine_global (1);
        
@@ -163,7 +163,7 @@ int main () {
       tria.execute_refinement ();
     };
   tria.refine_global (1);
-  
+*/
 
   cout << "Distributing dofs... "; 
   dof.distribute_dofs (fe);

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.