From ea979c657b12428aaf3fdc9ee9ba4347492a58e2 Mon Sep 17 00:00:00 2001 From: kanschat Date: Fri, 8 Oct 2010 20:48:45 +0000 Subject: [PATCH] automatically create quadrature if none is provided git-svn-id: https://svn.dealii.org/trunk@22295 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/mesh_worker_info.h | 287 ++++++++++------ deal.II/examples/step-39/step-39.cc | 40 +-- deal.II/lac/include/lac/matrix_block.h | 312 ++++++++++++++++-- 3 files changed, 491 insertions(+), 148 deletions(-) diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.h b/deal.II/deal.II/include/numerics/mesh_worker_info.h index ced843f5f9..f70e244997 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_info.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.h @@ -526,6 +526,10 @@ namespace MeshWorker *
  • It stores information on finite element vectors and whether * their data should be used to compute values or derivatives of * functions at quadrature points. + * + *
  • It makes educated guesses on quadrature rules and update + * flags, so that minimal code has to be written when default + * parameters are sufficient. * * * In order to allow for sufficient generality, a few steps have to be @@ -541,10 +545,16 @@ namespace MeshWorker * the selectors filled before and adds all the flags needed to get * the selection. Additional flags can be set with add_update_flags(). * - * Finally, we need to choose quadrature formulas. If you choose to - * use Gauss formulas only, use initialize_gauss_quadrature() with - * appropriate values. Otherwise, you can fill the variables - * #cell_quadrature, #boundary_quadrature and #face_quadrature directly. + * Finally, we need to choose quadrature formulas. In the simplest + * case, you might be happy with the default settings, which are + * n-point Gauss formulas. If only derivatives of the shape + * functions are used (#update_values is not set) n equals the + * highest polynomial degree in the FiniteElement, if #update_values + * is set, n is one higher than this degree. If you choose to + * use Gauss formulas of other size, use initialize_gauss_quadrature() + * with appropriate values. Otherwise, you can fill the variables + * #cell_quadrature, #boundary_quadrature and #face_quadrature + * directly. * * In order to save time, you can set the variables boundary_fluxes * and interior_fluxes of the base class to false, thus telling the @@ -571,17 +581,56 @@ namespace MeshWorker * Default constructor. */ IntegrationInfoBox (); - + + /** + * Initialize the + * IntegrationInfo objects + * contained. + * + * Before doing so, add update + * flags necessary to produce + * the data needed and also + * set uninitialized quadrature + * rules to Gauss formulas, + * which integrate polynomial + * bilinear forms exactly. + */ void initialize(const FiniteElement& el, const Mapping& mapping, const BlockInfo* block_info = 0); + /** + * Initialize the + * IntegrationInfo objects + * contained. + * + * Before doing so, add update + * flags necessary to produce + * the data needed and also + * set uninitialized quadrature + * rules to Gauss formulas, + * which integrate polynomial + * bilinear forms exactly. + */ template void initialize(const FiniteElement& el, const Mapping& mapping, const NamedData& data, const BlockInfo* block_info = 0); + /** + * Initialize the + * IntegrationInfo objects + * contained. + * + * Before doing so, add update + * flags necessary to produce + * the data needed and also + * set uninitialized quadrature + * rules to Gauss formulas, + * which integrate polynomial + * bilinear forms exactly. + */ template void initialize(const FiniteElement& el, const Mapping& mapping, @@ -637,17 +686,26 @@ namespace MeshWorker const bool face = true, const bool neighbor = true); - /** Assign n-point Gauss + /** + * Assign n-point Gauss * quadratures to each of the * quadrature rules. Here, a * size of zero points means * that no loop over these grid * entities should be * performed. + * + * If the parameter + * force is true, then + * all quadrature sets are + * filled with new quadrature + * ruels. If it is false, then + * only empty rules are changed. */ void initialize_gauss_quadrature(unsigned int n_cell_points, unsigned int n_boundary_points, - unsigned int n_face_points); + unsigned int n_face_points, + const bool force = true); /** * The memory used by this object. @@ -1303,42 +1361,105 @@ namespace MeshWorker } +//----------------------------------------------------------------------// + + template <> + inline + void + IntegrationInfoBox<1,1>::initialize_gauss_quadrature( + const unsigned int cp, + const unsigned int, + const unsigned int, + bool force) + { + if (force || cell_quadrature.size() == 0) + cell_quadrature = QGauss<1>(cp); + } + + + template <> + inline + void + IntegrationInfoBox<1,2>::initialize_gauss_quadrature( + const unsigned int cp, + const unsigned int, + const unsigned int, + bool force) + { + if (force || cell_quadrature.size() == 0) + cell_quadrature = QGauss<1>(cp); + } + + template inline void - IntegrationInfoBox:: - initialize(const FiniteElement& el, - const Mapping& mapping, - const BlockInfo* block_info) + IntegrationInfoBox::initialize_gauss_quadrature( + unsigned int cp, + unsigned int bp, + unsigned int fp, + bool force) { - initialize_update_flags(); - - cell.template initialize >(el, mapping, cell_quadrature, - cell_flags, block_info); - boundary.template initialize >(el, mapping, boundary_quadrature, - boundary_flags, block_info); - face.template initialize >(el, mapping, face_quadrature, - face_flags, block_info); - subface.template initialize >(el, mapping, face_quadrature, - face_flags, block_info); - neighbor.template initialize >(el, mapping, face_quadrature, - neighbor_flags, block_info); + if (force || cell_quadrature.size() == 0) + cell_quadrature = QGauss(cp); + if (force || boundary_quadrature.size() == 0) + boundary_quadrature = QGauss(bp); + if (force || face_quadrature.size() == 0) + face_quadrature = QGauss(fp); + } + + + template + inline + void + IntegrationInfoBox::add_update_flags_all (const UpdateFlags flags) + { + add_update_flags(flags, true, true, true, true); + } + + + template + inline + void + IntegrationInfoBox::add_update_flags_cell (const UpdateFlags flags) + { + add_update_flags(flags, true, false, false, false); } + + template + inline + void + IntegrationInfoBox::add_update_flags_boundary (const UpdateFlags flags) + { + add_update_flags(flags, false, true, false, false); + } + + template + inline + void + IntegrationInfoBox::add_update_flags_face (const UpdateFlags flags) + { + add_update_flags(flags, false, false, true, true); + } + + template <> inline void - IntegrationInfoBox<1,1>:: - initialize(const FiniteElement<1,1>& el, - const Mapping<1,1>& mapping, - const BlockInfo* block_info) + IntegrationInfoBox<1,1>::initialize( + const FiniteElement<1,1>& el, + const Mapping<1,1>& mapping, + const BlockInfo* block_info) { initialize_update_flags(); + initialize_gauss_quadrature( + (cell_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), 1, 1, false); const int dim = 1; const int sdim = 1; - + cell.initialize >(el, mapping, cell_quadrature, cell_flags, block_info); } @@ -1348,12 +1469,14 @@ namespace MeshWorker template <> inline void - IntegrationInfoBox<1,2>:: - initialize(const FiniteElement<1,2>& el, - const Mapping<1,2>& mapping, - const BlockInfo* block_info) + IntegrationInfoBox<1,2>::initialize( + const FiniteElement<1,2>& el, + const Mapping<1,2>& mapping, + const BlockInfo* block_info) { initialize_update_flags(); + initialize_gauss_quadrature( + (cell_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), 1, 1, false); const int dim = 1; const int sdim = 2; @@ -1363,7 +1486,32 @@ namespace MeshWorker } -//----------------------------------------------------------------------// + template + inline + void + IntegrationInfoBox::initialize( + const FiniteElement& el, + const Mapping& mapping, + const BlockInfo* block_info) + { + initialize_update_flags(); + initialize_gauss_quadrature( + (cell_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), + (boundary_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), + (face_flags & update_values) ? (el.tensor_degree()+1) : el.tensor_degree(), false); + + cell.template initialize >(el, mapping, cell_quadrature, + cell_flags, block_info); + boundary.template initialize >(el, mapping, boundary_quadrature, + boundary_flags, block_info); + face.template initialize >(el, mapping, face_quadrature, + face_flags, block_info); + subface.template initialize >(el, mapping, face_quadrature, + face_flags, block_info); + neighbor.template initialize >(el, mapping, face_quadrature, + neighbor_flags, block_info); + } + template template @@ -1441,79 +1589,6 @@ namespace MeshWorker {} - template - inline - void - IntegrationInfoBox::initialize_gauss_quadrature( - unsigned int cp, - unsigned int bp, - unsigned int fp) - { - cell_quadrature = QGauss(cp); - boundary_quadrature = QGauss(bp); - face_quadrature = QGauss(fp); - } - - - template <> - inline - void - IntegrationInfoBox<1,1>:: - initialize_gauss_quadrature(const unsigned int cp, - const unsigned int, - const unsigned int) - { - cell_quadrature = QGauss<1>(cp); - } - - - template <> - inline - void - IntegrationInfoBox<1,2>:: - initialize_gauss_quadrature(const unsigned int cp, - const unsigned int, - const unsigned int) - { - cell_quadrature = QGauss<1>(cp); - } - - - template - inline - void - IntegrationInfoBox::add_update_flags_all (const UpdateFlags flags) - { - add_update_flags(flags, true, true, true, true); - } - - - template - inline - void - IntegrationInfoBox::add_update_flags_cell (const UpdateFlags flags) - { - add_update_flags(flags, true, false, false, false); - } - - - template - inline - void - IntegrationInfoBox::add_update_flags_boundary (const UpdateFlags flags) - { - add_update_flags(flags, false, true, false, false); - } - - - template - inline - void - IntegrationInfoBox::add_update_flags_face (const UpdateFlags flags) - { - add_update_flags(flags, false, false, true, true); - } - } DEAL_II_NAMESPACE_CLOSE diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index 335bec6d78..e47b9deff3 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -570,27 +570,33 @@ Step39::assemble_matrix() // such that they always point to // the current cell. To this end, // we need to tell it first, where - // and what to compute, + // and what to compute. Since we + // are not doing anything fancy, we + // can rely on their standard + // choice for quadrature rules. + // + // Since their default update flags + // are minimal, we add what we need + // additionally, namely the values + // and gradients of shape functions + // on all objects (cells, boundary + // and interior faces). Afterwards, + // we are ready to initialize the + // container, which will create all + // necessary FEValuesBase objects + // for integration. MeshWorker::IntegrationInfoBox info_box; - // namely, which quadrature - // formulas to use and - const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; - info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); - // which values to update in these - // points. Update flags are - // initialized to some default - // values to be able to - // integrate. Here, we add what we - // need additionally, namely the - // values and gradients of shape - // functions on all objects (cells, - // boundary and interior faces). UpdateFlags update_flags = update_values | update_gradients; info_box.add_update_flags_all(update_flags); info_box.initialize(fe, mapping); // This is the object into which we - // integrate local data. + // integrate local data. It is + // filled by the local integration + // routines in MatrixIntegrator and + // then used by the assembler to + // distribute the information into + // the global matrix. MeshWorker::DoFInfo dof_info(dof_handler); // Finally, we need an object that @@ -628,8 +634,6 @@ void Step39::assemble_mg_matrix() { MeshWorker::IntegrationInfoBox info_box; - const unsigned int n_gauss_points = mg_dof_handler.get_fe().tensor_degree()+1; - info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points, n_gauss_points); UpdateFlags update_flags = update_values | update_gradients; info_box.add_update_flags_all(update_flags); info_box.initialize(fe, mapping); @@ -673,8 +677,6 @@ void Step39::assemble_right_hand_side() { MeshWorker::IntegrationInfoBox info_box; - const unsigned int n_gauss_points = dof_handler.get_fe().tensor_degree()+1; - info_box.initialize_gauss_quadrature(n_gauss_points, n_gauss_points+1, n_gauss_points); UpdateFlags update_flags = update_quadrature_points | update_values | update_gradients; info_box.add_update_flags_all(update_flags); info_box.initialize(fe, mapping); diff --git a/deal.II/lac/include/lac/matrix_block.h b/deal.II/lac/include/lac/matrix_block.h index 3560b7d993..1dba79fd2a 100644 --- a/deal.II/lac/include/lac/matrix_block.h +++ b/deal.II/lac/include/lac/matrix_block.h @@ -479,13 +479,35 @@ class MatrixBlockVector */ template class MGMatrixBlockVector - : private NamedData > > > + : public Subscriptor { public: /** * The type of object stored. */ typedef MGLevelObject > value_type; + /** + * Constructor, determining which + * matrices should be stored. + * + * If edge_matrices is + * true, then objects for edge + * matrices for discretizations + * with degrees of freedom on + * faces are allocated. + * + * If edge_flux_matrices + * is true, then objects for DG + * fluxes on the refinement edge + * are allocated. + */ + MGMatrixBlockVector(const bool edge_matrices = false, + const bool edge_flux_matrices = false); + + /** + * The number of blocks. + */ + unsigned int size () const; /** * Add a new matrix block at the @@ -504,8 +526,35 @@ class MGMatrixBlockVector * reinitializes each matrix in * the vector with the correct * pattern from the block system. + * + * This function reinitializes + * the level matrices. */ - void reinit(const MGLevelObject& sparsity); + void reinit_matrix(const MGLevelObject& sparsity); + /** + * For matrices using a + * SparsityPattern, this function + * reinitializes each matrix in + * the vector with the correct + * pattern from the block system. + * + * This function reinitializes + * the matrices for degrees of + * freedom on the refinement edge. + */ + void reinit_edge(const MGLevelObject& sparsity); + /** + * For matrices using a + * SparsityPattern, this function + * reinitializes each matrix in + * the vector with the correct + * pattern from the block system. + * + * This function reinitializes + * the flux matrices over the + * refinement edge. + */ + void reinit_edge_flux(const MGLevelObject& sparsity); /** * Clears the object. @@ -528,26 +577,88 @@ class MGMatrixBlockVector /** * Access a constant reference to - * the block at position i. + * the matrix block at position i. */ const value_type& block(unsigned int i) const; /** * Access a reference to - * the block at position i. + * the matrix block at position i. */ value_type& block(unsigned int i); + /** + * Access a constant reference to + * the edge matrix block at position i. + */ + const value_type& block_in(unsigned int i) const; + + /** + * Access a reference to + * the edge matrix block at position i. + */ + value_type& block_in(unsigned int i); + + /** + * Access a constant reference to + * the edge matrix block at position i. + */ + const value_type& block_out(unsigned int i) const; + + /** + * Access a reference to + * the edge matrix block at position i. + */ + value_type& block_out(unsigned int i); + + /** + * Access a constant reference to + * the edge flux matrix block at position i. + */ + const value_type& block_up(unsigned int i) const; + + /** + * Access a reference to + * the edge flux matrix block at position i. + */ + value_type& block_up(unsigned int i); + + /** + * Access a constant reference to + * the edge flux matrix block at position i. + */ + const value_type& block_down(unsigned int i) const; + + /** + * Access a reference to + * the edge flux matrix block at position i. + */ + value_type& block_down(unsigned int i); + /** * The memory used by this object. */ unsigned int memory_consumption () const; - - NamedData >::subscribe; - NamedData >::unsubscribe; - NamedData >::size; - private: + /// Clear one of the matrix objects + void clear_object(NamedData > >&); + + /// Flag for storing #matrices_in and #matrices_out + const bool edge_matrices; + + /// Flag for storing #flux_matrices_up and #flux_matrices_down + const bool edge_flux_matrices; + + /// The level matrices + NamedData > > matrices; + /// The matrix from the interior of a level to the refinement edge + NamedData > > matrices_in; + /// The matrix from the refinement edge to the interior of a level + NamedData > > matrices_out; + /// The DG flux from a level to the lower level + NamedData > > flux_matrices_down; + /// The DG flux from the lower level to a level + NamedData > > flux_matrices_up; }; @@ -871,18 +982,46 @@ MatrixBlockVector::matrix(unsigned int i) //----------------------------------------------------------------------// +template +inline +MGMatrixBlockVector::MGMatrixBlockVector( + const bool e, const bool f) + : + edge_matrices(e), + edge_flux_matrices(f) +{} + + +template +inline +unsigned int +MGMatrixBlockVector::size () const +{ + return matrices.size(); +} + + template inline void MGMatrixBlockVector::add( unsigned int row, unsigned int column, const std::string& name) { - std_cxx1x::shared_ptr > > - p(new value_type(0, 1)); - (*p)[0].row = row; - (*p)[0].column = column; + MGLevelObject > p(0, 1); + p[0].row = row; + p[0].column = column; - NamedData >::add(p, name); + matrices.add(p, name); + if (edge_matrices) + { + matrices_in.add(p, name); + matrices_out.add(p, name); + } + if (edge_flux_matrices) + { + flux_matrices_up.add(p, name); + flux_matrices_down.add(p, name); + } } @@ -890,7 +1029,7 @@ template inline const MGLevelObject >& MGMatrixBlockVector::block(unsigned int i) const { - return *this->read(i); + return matrices.read(i); } @@ -898,13 +1037,77 @@ template inline MGLevelObject >& MGMatrixBlockVector::block(unsigned int i) { - return *(*this)(i); + return matrices(i); +} + + +template +inline const MGLevelObject >& +MGMatrixBlockVector::block_in(unsigned int i) const +{ + return matrices_in.read(i); +} + + +template +inline MGLevelObject >& +MGMatrixBlockVector::block_in(unsigned int i) +{ + return matrices_in(i); +} + + +template +inline const MGLevelObject >& +MGMatrixBlockVector::block_out(unsigned int i) const +{ + return matrices_out.read(i); +} + + +template +inline MGLevelObject >& +MGMatrixBlockVector::block_out(unsigned int i) +{ + return matrices_out(i); +} + + +template +inline const MGLevelObject >& +MGMatrixBlockVector::block_up(unsigned int i) const +{ + return flux_matrices_up.read(i); +} + + +template +inline MGLevelObject >& +MGMatrixBlockVector::block_up(unsigned int i) +{ + return flux_matrices_up(i); +} + + +template +inline const MGLevelObject >& +MGMatrixBlockVector::block_down(unsigned int i) const +{ + return flux_matrices_down.read(i); +} + + +template +inline MGLevelObject >& +MGMatrixBlockVector::block_down(unsigned int i) +{ + return flux_matrices_down(i); } template inline void -MGMatrixBlockVector::reinit(const MGLevelObject& sparsity) +MGMatrixBlockVector::reinit_matrix(const MGLevelObject& sparsity) { for (unsigned int i=0;isize();++i) { @@ -923,6 +1126,70 @@ MGMatrixBlockVector::reinit(const MGLevelObject& s } +template +inline void +MGMatrixBlockVector::reinit_edge(const MGLevelObject& sparsity) +{ + for (unsigned int i=0;isize();++i) + { + MGLevelObject >& o = block(i); + const unsigned int row = o[o.min_level()].row; + const unsigned int col = o[o.min_level()].column; + + block_in(i).resize(sparsity.min_level(), sparsity.max_level()); + block_out(i).resize(sparsity.min_level(), sparsity.max_level()); + for (unsigned int level = o.min_level();level <= o.max_level();++level) + { + block_in(i)[level].row = row; + block_in(i)[level].column = col; + internal::reinit(block_in(i)[level], sparsity[level]); + block_out(i)[level].row = row; + block_out(i)[level].column = col; + internal::reinit(block_out(i)[level], sparsity[level]); + } + } +} + + +template +inline void +MGMatrixBlockVector::reinit_edge_flux(const MGLevelObject& sparsity) +{ + for (unsigned int i=0;isize();++i) + { + MGLevelObject >& o = block(i); + const unsigned int row = o[o.min_level()].row; + const unsigned int col = o[o.min_level()].column; + + block_up(i).resize(sparsity.min_level(), sparsity.max_level()); + block_down(i).resize(sparsity.min_level(), sparsity.max_level()); + for (unsigned int level = o.min_level();level <= o.max_level();++level) + { + block_up(i)[level].row = row; + block_up(i)[level].column = col; + internal::reinit(block_up(i)[level], sparsity[level]); + block_down(i)[level].row = row; + block_down(i)[level].column = col; + internal::reinit(block_down(i)[level], sparsity[level]); + } + + } +} + + +template +inline void +MGMatrixBlockVector::clear_object(NamedData > >& mo) +{ + for (unsigned int i=0;i >& o = mo(i); + for (unsigned int level = o.min_level();level <= o.max_level();++level) + o[level].matrix.clear(); + } +} + + template inline void MGMatrixBlockVector::clear(bool really_clean) @@ -933,12 +1200,11 @@ MGMatrixBlockVector::clear(bool really_clean) } else { - for (unsigned int i=0;isize();++i) - { - MGLevelObject >& o = block(i); - for (unsigned int level = o.min_level();level <= o.max_level();++level) - o[level].matrix.clear(); - } + clear_object(matrices); + clear_object(matrices_in); + clear_object(matrices_out); + clear_object(flux_matrices_up); + clear_object(flux_matrices_down); } } -- 2.39.5