From b7ba5f591b5bc1c1263e7efb6d8b74b0a6b71c44 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 24 Apr 1998 08:12:58 +0000 Subject: [PATCH] Change the way we specify which fields in FEValues need to be updated. git-svn-id: https://svn.dealii.org/trunk@190 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/Attic/examples/Makefile | 16 ++- .../Attic/examples/convergence/convergence.cc | 7 +- .../deal.II/Attic/examples/poisson/problem.cc | 5 +- deal.II/deal.II/include/fe/fe.h | 117 +++++++++--------- deal.II/deal.II/include/numerics/assembler.h | 4 +- deal.II/deal.II/include/numerics/base.h | 10 +- deal.II/deal.II/source/fe/fe.cc | 40 +++--- deal.II/deal.II/source/numerics/assembler.cc | 2 +- deal.II/deal.II/source/numerics/base.cc | 19 ++- tests/big-tests/Makefile | 16 ++- tests/big-tests/convergence/convergence.cc | 7 +- tests/big-tests/poisson/problem.cc | 5 +- 12 files changed, 124 insertions(+), 124 deletions(-) diff --git a/deal.II/deal.II/Attic/examples/Makefile b/deal.II/deal.II/Attic/examples/Makefile index 9d1e1501ce..fb4b0f3c29 100644 --- a/deal.II/deal.II/Attic/examples/Makefile +++ b/deal.II/deal.II/Attic/examples/Makefile @@ -25,7 +25,8 @@ ifeq ($(shell uname),SunOS) CXX = /usr/local/gcc/gcc-2.8.1/bin/c++ endif -cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc convergence/convergence.cc +cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc \ + convergence/convergence.cc error-estimation/error-estimation o-files = $(cc-files:.cc=.o) h-files = $(wildcard ../include/*.h) @@ -50,8 +51,12 @@ convergence/convergence: convergence/convergence.o $(LIBFILES.g) @echo ================= Linking $@ @$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o +error-estimation/error-estimation: error-estimation/error-estimation.o $(LIBFILES.g) + @echo ================= Linking $@ + @$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o -run: run_grid_test run_dof_test run_poisson_test run_convergence_test + +run: run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test run_grid_test: cd grid ; grid_test 4 ; mv grid.1 grid.4 @@ -71,14 +76,19 @@ run_convergence_test: cd convergence ; convergence cd convergence ; gnuplot make_ps +run_error_test: + cd error-estimation ; error-estimation + cd error-estimation ; gnuplot make_ps + clean: cd grid ; rm -f grid.[1234] *.eps *.o *~ grid_test cd dof ; rm -f grid.* sparsity.* *.o *~ dof_test cd poisson ; make clean cd convergence ; rm -f *.o convergence *~ *.eps gnuplot.* + cd error-estimation ; rm -f *.o error-estimation *~ *.eps gnuplot.* -.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test clean +.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test clean #Rule to generate the dependency file. This file is diff --git a/deal.II/deal.II/Attic/examples/convergence/convergence.cc b/deal.II/deal.II/Attic/examples/convergence/convergence.cc index cff66569ed..ecea08cec3 100644 --- a/deal.II/deal.II/Attic/examples/convergence/convergence.cc +++ b/deal.II/deal.II/Attic/examples/convergence/convergence.cc @@ -251,9 +251,8 @@ void PoissonProblem::run (const unsigned int level) { n_dofs.push_back (dof->n_dofs()); cout << " Assembling matrices..." << endl; - FEValues::UpdateStruct update_flags; - update_flags.q_points = update_flags.gradients = true; - update_flags.jacobians = update_flags.JxW_values = true; + UpdateFields update_flags = UpdateFields(update_q_points | update_gradients | + update_jacobians | update_JxW_values); ProblemBase::DirichletBC dirichlet_bc; dirichlet_bc[0] = boundary_values; @@ -372,7 +371,7 @@ void PoissonProblem::print_history () const { int main () { PoissonProblem<2> problem; - for (unsigned int level=1; level<10; ++level) + for (unsigned int level=1; level<9; ++level) problem.run (level); cout << endl << "Printing convergence history to ..." << endl; diff --git a/deal.II/deal.II/Attic/examples/poisson/problem.cc b/deal.II/deal.II/Attic/examples/poisson/problem.cc index 8a7e818576..553a0c36cd 100644 --- a/deal.II/deal.II/Attic/examples/poisson/problem.cc +++ b/deal.II/deal.II/Attic/examples/poisson/problem.cc @@ -432,10 +432,7 @@ void PoissonProblem::run (ParameterHandler &prm) { cout << dof->n_dofs() << " degrees of freedom." << endl; cout << " Assembling matrices..." << endl; - FEValues::UpdateStruct update_flags; - update_flags.q_points = update_flags.gradients = true; - update_flags.jacobians = update_flags.JxW_values = true; - + UpdateFields update_flags = UpdateFields(update_gradients | update_JxW_values); ProblemBase::DirichletBC dirichlet_bc; dirichlet_bc[0] = boundary_values; assemble (equation, quadrature, fe, update_flags, dirichlet_bc); diff --git a/deal.II/deal.II/include/fe/fe.h b/deal.II/deal.II/include/fe/fe.h index 73fcd32aaf..482e694caa 100644 --- a/deal.II/deal.II/include/fe/fe.h +++ b/deal.II/deal.II/include/fe/fe.h @@ -16,6 +16,54 @@ template class FiniteElement; template class Quadrature; +/** + Provide a set of flags which tells the #FEValues<>::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 variable of this type has to be passed to the constructor of the + #FEValues# object. You can select more than one flag by concatenation + using the #|# (bitwise #or#) operator. + */ +enum UpdateFields { + /** + * Default: update nothing. + */ + update_default = 0, + /** + * Compute quadrature points in real + * space (not on unit cell). + */ + update_q_points = 1, + /** + * Transform gradients on unit cell to + * gradients on real cell. + */ + update_gradients = 2, + /** + * Compute jacobian matrices of the + * transform between unit and real cell + * in the evaluation points. + */ + update_jacobians = 4, + /** + * Compute the JxW values (Jacobian + * determinant at the quadrature point + * times the weight of this point). + */ + update_JxW_values = 8, + /** + * Compute the points on the real cell + * on which the ansatz functions are + * located. + */ + update_ansatz_points = 16 +}; + + + + /** Represent a finite element evaluated with a specific quadrature rule. This class is an optimization which avoids evaluating the shape functions @@ -28,7 +76,7 @@ template class Quadrature; The unit cell is defined to be the tensor product of the interval $[0,1]$ in the present number of dimensions. In part of the literature, the convention - is used that the unit cell be the tensor product of the intervale $[-1,1]$, + is used that the unit cell be the tensor product of the interval $[-1,1]$, which is to distinguished properly. Objects of this class store a multitude of different values needed to @@ -41,7 +89,9 @@ template class Quadrature; The Jacobian matrix is defined to be $$ J_{ij} = {d\xi_i \over dx_j} $$ - which is the form needed to compute the gradient on the real cell from + where the $\xi_i$ are the coordinates on the unit cell and the $x_i$ are + the coordinates on the real cell. + This is the form needed to compute the gradient on the real cell from the gradient on the unit cell. If we want to transform the area element $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 @@ -50,7 +100,7 @@ template class Quadrature; 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 + object knows about which fields are needed by the #UpdateFields# 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. @@ -58,55 +108,6 @@ template class Quadrature; 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 q_points; - /** - * Transform gradients on unit cell to - * gradients on real cell. - */ - bool gradients; - /** - * Compute jacobian matrices of the - * transform between unit and real cell - * in the evaluation points. - */ - bool jacobians; - /** - * Compute the JxW values (Jacobian - * determinant at the quadrature point - * times the weight of this point). - */ - bool JxW_values; - /** - * Compute the points on the real cell - * on which the ansatz functions are - * located. - */ - bool ansatz_points; - }; @@ -129,7 +130,7 @@ class FEValues { */ FEValues (const FiniteElement &, const Quadrature &, - const UpdateStruct &); + const UpdateFields); /** * Return the value of the #i#th shape @@ -339,7 +340,7 @@ class FEValues { * Store which fields are to be updated by * the reinit function. */ - UpdateStruct update_flags; + UpdateFields update_flags; }; @@ -907,7 +908,7 @@ template inline const vector > > & FEValues::get_shape_grads () const { - Assert (update_flags.gradients, ExcAccessToUninitializedField()); + Assert (update_flags | update_gradients, ExcAccessToUninitializedField()); return shape_gradients; }; @@ -917,7 +918,7 @@ template inline const vector > & FEValues::get_quadrature_points () const { - Assert (update_flags.q_points, ExcAccessToUninitializedField()); + Assert (update_flags | update_q_points, ExcAccessToUninitializedField()); return quadrature_points; }; @@ -927,7 +928,7 @@ template inline const vector > & FEValues::get_ansatz_points () const { - Assert (update_flags.ansatz_points, ExcAccessToUninitializedField()); + Assert (update_flags | update_ansatz_points, ExcAccessToUninitializedField()); return ansatz_points; }; @@ -937,7 +938,7 @@ template inline const vector & FEValues::get_JxW_values () const { - Assert (update_flags.JxW_values, ExcAccessToUninitializedField()); + Assert (update_flags | update_JxW_values, ExcAccessToUninitializedField()); return JxW_values; }; diff --git a/deal.II/deal.II/include/numerics/assembler.h b/deal.II/deal.II/include/numerics/assembler.h index 0227875c11..f8c5f78888 100644 --- a/deal.II/deal.II/include/numerics/assembler.h +++ b/deal.II/deal.II/include/numerics/assembler.h @@ -106,7 +106,7 @@ struct AssemblerData { dVector &rhs_vector, const Quadrature &quadrature, const FiniteElement &fe, - const FEValues::UpdateStruct &update_flags); + const UpdateFields &update_flags); /** * Pointer to the dof handler object @@ -151,7 +151,7 @@ struct AssemblerData { * FEValues object need to be reinitialized * on each cell. */ - const FEValues::UpdateStruct update_flags; + const UpdateFields update_flags; }; diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 0e04b3cef3..0b943cfd88 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -290,11 +290,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 FEValues::UpdateStruct &update_flags, - const DirichletBC &dirichlet_bc = DirichletBC()); + virtual void assemble (const Equation &equation, + const Quadrature &q, + const FiniteElement &fe, + const UpdateFields &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 7425da5514..b28f842584 100644 --- a/deal.II/deal.II/source/fe/fe.cc +++ b/deal.II/deal.II/source/fe/fe.cc @@ -7,16 +7,6 @@ -template -FEValues::UpdateStruct::UpdateStruct () : - q_points(false), - gradients(false), - jacobians(false), - JxW_values(false), - ansatz_points(false) {}; - - - /*------------------------------- FEValues -------------------------------*/ @@ -25,7 +15,7 @@ FEValues::UpdateStruct::UpdateStruct () : template FEValues::FEValues (const FiniteElement &fe, const Quadrature &quadrature, - const UpdateStruct &update_flags) : + const UpdateFields update_flags) : n_quadrature_points(quadrature.n_quadrature_points), total_dofs(fe.total_dofs), shape_values(fe.total_dofs, quadrature.n_quadrature_points), @@ -76,7 +66,7 @@ FEValues::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.gradients, ExcAccessToUninitializedField()); + Assert (update_flags | update_gradients, ExcAccessToUninitializedField()); return shape_gradients[i][j]; }; @@ -86,7 +76,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 const Point & FEValues::ansatz_point (const unsigned int i) const { Assert (i & FEValues::ansatz_point (const unsigned int i) const { template double FEValues::JxW (const unsigned int i) const { Assert (i::reinit (const typename Triangulation::cell_iterator &ce const FiniteElement &fe) { // fill jacobi matrices and real // quadrature points - if (update_flags.jacobians || update_flags.q_points) + if ((update_flags | update_jacobians) || + (update_flags | update_q_points)) fe.fill_fe_values (cell, unit_quadrature_points, jacobi_matrices, - update_flags.jacobians, + update_flags | update_jacobians, ansatz_points, - update_flags.ansatz_points, + update_flags | update_ansatz_points, quadrature_points, - update_flags.q_points); + update_flags | update_q_points); // compute gradients on real element if // requested - if (update_flags.gradients) + if (update_flags | update_gradients) { - Assert (update_flags.jacobians, ExcCannotInitializeField()); + Assert (update_flags | update_jacobians, ExcCannotInitializeField()); for (unsigned int i=0; i::reinit (const typename Triangulation::cell_iterator &ce // refer to the general doc for // why we take the inverse of the // determinant - if (update_flags.JxW_values) + if (update_flags | update_JxW_values) { - Assert (update_flags.jacobians, + Assert (update_flags | update_jacobians, ExcCannotInitializeField()); for (unsigned int i=0; i::face_ansatz_points (const typename Triangulation<2>::face /*------------------------------- 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/numerics/assembler.cc b/deal.II/deal.II/source/numerics/assembler.cc index 66411ddf05..cf47c72386 100644 --- a/deal.II/deal.II/source/numerics/assembler.cc +++ b/deal.II/deal.II/source/numerics/assembler.cc @@ -51,7 +51,7 @@ AssemblerData::AssemblerData (const DoFHandler &dof, dVector &rhs_vector, const Quadrature &quadrature, const FiniteElement &fe, - const FEValues::UpdateStruct &update_flags) : + const UpdateFields &update_flags) : dof(dof), assemble_matrix(assemble_matrix), assemble_rhs(assemble_rhs), diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index a905423994..e5b3cd1875 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -73,11 +73,11 @@ ProblemBase::~ProblemBase () {}; template -void ProblemBase::assemble (const Equation &equation, - const Quadrature &quadrature, - const FiniteElement &fe, - const FEValues::UpdateStruct &update_flags, - const DirichletBC &dirichlet_bc) { +void ProblemBase::assemble (const Equation &equation, + const Quadrature &quadrature, + const FiniteElement &fe, + const UpdateFields &update_flags, + const DirichletBC &dirichlet_bc) { Assert ((tria!=0) && (dof_handler!=0), ExcNoTriaSelected()); system_sparsity.reinit (dof_handler->n_dofs(), @@ -164,12 +164,11 @@ void ProblemBase::integrate_difference (const Function &exact_sol difference.reinit (tria->n_active_cells()); - FEValues::UpdateStruct update_flags; - update_flags.q_points = true; - update_flags.jacobians = true; - update_flags.JxW_values = true; + UpdateFields update_flags = UpdateFields (update_q_points | + update_jacobians | + update_JxW_values); if ((norm==H1_seminorm) || (norm==H1_norm)) - update_flags.gradients = true; + update_flags = UpdateFields (update_flags | update_gradients); FEValues fe_values(fe, q, update_flags); // loop over all cells diff --git a/tests/big-tests/Makefile b/tests/big-tests/Makefile index 9d1e1501ce..fb4b0f3c29 100644 --- a/tests/big-tests/Makefile +++ b/tests/big-tests/Makefile @@ -25,7 +25,8 @@ ifeq ($(shell uname),SunOS) CXX = /usr/local/gcc/gcc-2.8.1/bin/c++ endif -cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc convergence/convergence.cc +cc-files = grid/grid_test.cc dof/dof_test.cc poisson/poisson.cc \ + convergence/convergence.cc error-estimation/error-estimation o-files = $(cc-files:.cc=.o) h-files = $(wildcard ../include/*.h) @@ -50,8 +51,12 @@ convergence/convergence: convergence/convergence.o $(LIBFILES.g) @echo ================= Linking $@ @$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o +error-estimation/error-estimation: error-estimation/error-estimation.o $(LIBFILES.g) + @echo ================= Linking $@ + @$(CXX) $(CXXFLAGS.g) -o $@ $< $(LIBS) ../../mia/control.o -run: run_grid_test run_dof_test run_poisson_test run_convergence_test + +run: run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test run_grid_test: cd grid ; grid_test 4 ; mv grid.1 grid.4 @@ -71,14 +76,19 @@ run_convergence_test: cd convergence ; convergence cd convergence ; gnuplot make_ps +run_error_test: + cd error-estimation ; error-estimation + cd error-estimation ; gnuplot make_ps + clean: cd grid ; rm -f grid.[1234] *.eps *.o *~ grid_test cd dof ; rm -f grid.* sparsity.* *.o *~ dof_test cd poisson ; make clean cd convergence ; rm -f *.o convergence *~ *.eps gnuplot.* + cd error-estimation ; rm -f *.o error-estimation *~ *.eps gnuplot.* -.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test clean +.PHONY: run run_grid_test run_dof_test run_poisson_test run_convergence_test run_error_test clean #Rule to generate the dependency file. This file is diff --git a/tests/big-tests/convergence/convergence.cc b/tests/big-tests/convergence/convergence.cc index cff66569ed..ecea08cec3 100644 --- a/tests/big-tests/convergence/convergence.cc +++ b/tests/big-tests/convergence/convergence.cc @@ -251,9 +251,8 @@ void PoissonProblem::run (const unsigned int level) { n_dofs.push_back (dof->n_dofs()); cout << " Assembling matrices..." << endl; - FEValues::UpdateStruct update_flags; - update_flags.q_points = update_flags.gradients = true; - update_flags.jacobians = update_flags.JxW_values = true; + UpdateFields update_flags = UpdateFields(update_q_points | update_gradients | + update_jacobians | update_JxW_values); ProblemBase::DirichletBC dirichlet_bc; dirichlet_bc[0] = boundary_values; @@ -372,7 +371,7 @@ void PoissonProblem::print_history () const { int main () { PoissonProblem<2> problem; - for (unsigned int level=1; level<10; ++level) + for (unsigned int level=1; level<9; ++level) problem.run (level); cout << endl << "Printing convergence history to ..." << endl; diff --git a/tests/big-tests/poisson/problem.cc b/tests/big-tests/poisson/problem.cc index 8a7e818576..553a0c36cd 100644 --- a/tests/big-tests/poisson/problem.cc +++ b/tests/big-tests/poisson/problem.cc @@ -432,10 +432,7 @@ void PoissonProblem::run (ParameterHandler &prm) { cout << dof->n_dofs() << " degrees of freedom." << endl; cout << " Assembling matrices..." << endl; - FEValues::UpdateStruct update_flags; - update_flags.q_points = update_flags.gradients = true; - update_flags.jacobians = update_flags.JxW_values = true; - + UpdateFields update_flags = UpdateFields(update_gradients | update_JxW_values); ProblemBase::DirichletBC dirichlet_bc; dirichlet_bc[0] = boundary_values; assemble (equation, quadrature, fe, update_flags, dirichlet_bc); -- 2.39.5