From f11068359c42a0e270dea06237a7c850f0095cc9 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 1 Apr 1998 12:37:54 +0000 Subject: [PATCH] More numerics and more efficient version of FEValues::reinit git-svn-id: https://svn.dealii.org/trunk@108 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/Attic/examples/poisson/poisson.cc | 14 +-- deal.II/deal.II/include/fe/fe.h | 100 ++++++++++++++++- deal.II/deal.II/include/fe/fe_lib.lagrange.h | 12 +- deal.II/deal.II/include/numerics/assembler.h | 9 +- deal.II/deal.II/include/numerics/base.h | 9 +- deal.II/deal.II/source/fe/fe.cc | 103 ++++++++++++----- deal.II/deal.II/source/fe/fe_lib.linear.cc | 104 +++++++++++------- deal.II/deal.II/source/numerics/assembler.cc | 9 +- deal.II/deal.II/source/numerics/base.cc | 25 +++-- tests/big-tests/poisson/poisson.cc | 14 +-- 10 files changed, 290 insertions(+), 109 deletions(-) diff --git a/deal.II/deal.II/Attic/examples/poisson/poisson.cc b/deal.II/deal.II/Attic/examples/poisson/poisson.cc index a1f6444408..03fa5a31ab 100644 --- a/deal.II/deal.II/Attic/examples/poisson/poisson.cc +++ b/deal.II/deal.II/Attic/examples/poisson/poisson.cc @@ -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); diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index f74d3e8414..f3ffc7d85c 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -41,10 +41,70 @@ template 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 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 &, - const Quadrature &); + const Quadrature &, + 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 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::cell_iterator &cell, const vector > &unit_points, vector &jacobians, - vector > &points) const; + const bool compute_jacobians, + vector > &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 > &unit_points, vector &jacobians, - vector > &points) const; + const bool compute_jacobians, + vector > &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 > &unit_points, vector &jacobians, - vector > &points) const; + const bool compute_jacobians, + vector > &points, + const bool compute_points) const; }; @@ -691,6 +780,7 @@ template inline const vector > > & FEValues::get_shape_grads () const { + Assert (update_flags.update_gradients, ExcAccessToUninitializedField()); return shape_gradients; }; @@ -700,6 +790,7 @@ template inline const vector > & FEValues::get_quadrature_points () const { + Assert (update_flags.update_q_points, ExcAccessToUninitializedField()); return quadrature_points; }; @@ -709,6 +800,7 @@ template inline const vector & FEValues::get_JxW_values () const { + Assert (update_flags.update_JxW_values, ExcAccessToUninitializedField()); return JxW_values; }; diff --git a/deal.II/deal.II/include/fe/fe_lib.lagrange.h b/deal.II/deal.II/include/fe/fe_lib.lagrange.h index e14f219e40..279f872740 100644 --- a/deal.II/deal.II/include/fe/fe_lib.lagrange.h +++ b/deal.II/deal.II/include/fe/fe_lib.lagrange.h @@ -70,7 +70,9 @@ class FELinear : public FiniteElement { virtual void fill_fe_values (const Triangulation::cell_iterator &cell, const vector > &unit_points, vector &jacobians, - vector > &points) const; + const bool compute_jacobians, + vector > &points, + const bool compute_points) const; }; @@ -129,7 +131,9 @@ class FEQuadratic : public FiniteElement { virtual void fill_fe_values (const Triangulation::cell_iterator &cell, const vector > &unit_points, vector &jacobians, - vector > &points) const; + const bool compute_jacobians, + vector > &points, + const bool compute_points) const; }; @@ -188,7 +192,9 @@ class FECubic : public FiniteElement { virtual void fill_fe_values (const Triangulation::cell_iterator &cell, const vector > &unit_points, vector &jacobians, - vector > &points) const; + const bool compute_jacobians, + vector > &points, + const bool compute_points) const; }; diff --git a/deal.II/deal.II/include/numerics/assembler.h b/deal.II/deal.II/include/numerics/assembler.h index bd2192629f..44af7879eb 100644 --- a/deal.II/deal.II/include/numerics/assembler.h +++ b/deal.II/deal.II/include/numerics/assembler.h @@ -105,7 +105,8 @@ struct AssemblerData { const bool assemble_rhs, ProblemBase &problem, const Quadrature &quadrature, - const FiniteElement &fe); + const FiniteElement &fe, + const FEValues::UpdateStruct &update_flags); /** * Pointer to the dof handler object @@ -135,6 +136,12 @@ struct AssemblerData { * object. */ const FiniteElement &fe; + /** + * Store which of the fields of the + * FEValues object need to be reinitialized + * on each cell. + */ + const FEValues::UpdateStruct update_flags; }; diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 572239610d..cfc2d4297d 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -248,10 +248,11 @@ class ProblemBase { * For what exactly happens here, refer to * the general doc of this class. */ - virtual void assemble (const Equation &equation, - const Quadrature &q, - const FiniteElement &fe, - const DirichletBC &dirichlet_bc = DirichletBC()); + virtual void assemble (const Equation &equation, + const Quadrature &q, + const FiniteElement &fe, + const FEValues::UpdateStruct &update_flags, + const DirichletBC &dirichlet_bc = DirichletBC()); /** * Solve the system of equations. diff --git a/deal.II/deal.II/source/fe/fe.cc b/deal.II/deal.II/source/fe/fe.cc index 7c591c7c4e..8e0895ebc0 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -7,12 +7,25 @@ +template +FEValues::UpdateStruct::UpdateStruct () : + update_q_points(false), + update_gradients(false), + update_jacobians(false), + update_JxW_values(false), + update_ansatz_points(false) {}; + + + + + /*------------------------------- FEValues -------------------------------*/ template FEValues::FEValues (const FiniteElement &fe, - const Quadrature &quadrature) : + const Quadrature &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::FEValues (const FiniteElement &fe, quadrature_points(quadrature.n_quadrature_points, Point()), unit_quadrature_points(quadrature.n_quadrature_points, Point()), jacobi_matrices (quadrature.n_quadrature_points, - dFMatrix(dim,dim)) + dFMatrix(dim,dim)), + update_flags (update_flags) { for (unsigned int i=0; i::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::shape_grad (const unsigned int i, template const Point & FEValues::quadrature_point (const unsigned int i) const { Assert (i & FEValues::quadrature_point (const unsigned int i) const template double FEValues::JxW (const unsigned int i) const { Assert (i void FEValues::reinit (const typename Triangulation::cell_iterator &cell, const FiniteElement &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 void FiniteElementBase::fill_fe_values (const typename Triangulation::cell_iterator &, const vector > &, vector &, - vector > &) const { + const bool, + vector > &, + 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 > &unit_points, vector &jacobians, - vector > &points) const { + const bool compute_jacobians, + vector > &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; ivertex(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 > &, vector &, - vector > &) const { + const bool, + vector > &, + 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>; diff --git a/deal.II/deal.II/source/fe/fe_lib.linear.cc b/deal.II/deal.II/source/fe/fe_lib.linear.cc index d0a167698b..a9c6fd16a2 100644 --- a/deal.II/deal.II/source/fe/fe_lib.linear.cc +++ b/deal.II/deal.II/source/fe/fe_lib.linear.cc @@ -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 > &unit_points, vector &jacobians, - vector > &points) const { + const bool compute_jacobians, + vector > &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 > &unit_points, vector &jacobians, - vector > &points) const { + const bool compute_jacobians, + vector > &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 vertices[n_vertices]; for (unsigned int l=0; lvertex(l); + + if (compute_points) + { + // initialize points to zero + for (unsigned int i=0; i (); + + // 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 (); - - // 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::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 gradient - = FELinear::shape_grad (s, unit_points[l]); - for (unsigned int i=0; i gradient + = FELinear::shape_grad (s, unit_points[l]); + for (unsigned int i=0; i::FEQuadratic () : void FEQuadratic<1>::fill_fe_values (const Triangulation<1>::cell_iterator &cell, const vector > &unit_points, vector &jacobians, - vector > &points) const { + const bool compute_jacobians, + vector > &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::shape_grad (const unsigned int i, void FEQuadratic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, const vector > &, vector &, - vector > &) const { + const bool, + vector > &, + 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 > &unit_points, vector &jacobians, - vector > &points) const { + const bool compute_jacobians, + vector > &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::shape_grad (const unsigned int i, void FECubic<2>::fill_fe_values (const Triangulation<2>::cell_iterator &, const vector > &, vector &, - vector > &) const { + const bool, + vector > &, + const bool) const { Assert (false, typename FiniteElementBase<2>::ExcNotImplemented()); }; diff --git a/deal.II/deal.II/source/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc index 218a66cddb..591e2de10d 100644 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ b/deal.II/deal.II/source/numerics/assembler.cc @@ -50,13 +50,15 @@ AssemblerData::AssemblerData (const DoFHandler &dof, const bool assemble_rhs, ProblemBase &problem, const Quadrature &quadrature, - const FiniteElement &fe) : + const FiniteElement &fe, + const FEValues::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::Assembler (Triangulation *tria, problem (((AssemblerData*)local_data)->problem), fe(((AssemblerData*)local_data)->fe), fe_values (((AssemblerData*)local_data)->fe, - ((AssemblerData*)local_data)->quadrature) + ((AssemblerData*)local_data)->quadrature, + ((AssemblerData*)local_data)->update_flags) { Assert ((unsigned int)problem.system_matrix.m() == dof_handler->n_dofs(), ExcInvalidData()); diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 84c946b708..64bc08bfc7 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -45,10 +45,11 @@ ProblemBase::~ProblemBase () {}; template -void ProblemBase::assemble (const Equation &equation, - const Quadrature &quadrature, - const FiniteElement &fe, - const DirichletBC &dirichlet_bc) { +void ProblemBase::assemble (const Equation &equation, + const Quadrature &quadrature, + const FiniteElement &fe, + const FEValues::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::assemble (const Equation &equation, // reinit solution vector, preset // with zeroes. solution.reinit (dof_handler->n_dofs()); - + // create assembler AssemblerData data (*dof_handler, true, true, //assemble matrix and rhs *this, quadrature, - fe); + fe, + update_flags); TriaActiveIterator > 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::integrate_difference (const Function &exact_sol difference.erase (difference.begin(), difference.end()); difference.reserve (tria->n_cells()); - FEValues fe_values(fe, q); + FEValues::UpdateStruct update_flags; + update_flags.update_q_points = true; + update_flags.update_JxW_values = true; + FEValues fe_values(fe, q, update_flags); // loop over all cells // (we need an iterator on the triangulation for @@ -187,9 +191,10 @@ void ProblemBase::integrate_difference (const Function &exact_sol psi); // then subtract finite element // solution - for (unsigned int j=0; j &JxW_values = fe_values.get_JxW_values(); + for (unsigned int j=0; j(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); -- 2.39.5