From d428c88234db8d5ee55cbbc767ae04356b533128 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Wed, 10 Feb 2010 03:06:44 +0000 Subject: [PATCH] change function calls and data structures in mesh worker to conform to work stream git-svn-id: https://svn.dealii.org/trunk@20536 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/numerics/mesh_worker_assembler.h | 90 ++--- .../include/numerics/mesh_worker_info.h | 331 +++++++++++------- .../numerics/mesh_worker_info.templates.h | 24 +- .../include/numerics/mesh_worker_loop.h | 71 ++-- .../source/numerics/mesh_worker_info.cc | 8 +- deal.II/examples/step-38/step-38.cc | 39 ++- deal.II/examples/step-39/step-39.cc | 121 ++++--- 7 files changed, 396 insertions(+), 288 deletions(-) diff --git a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h index f330e125c7..8eecafe100 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_assembler.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_assembler.h @@ -995,7 +995,7 @@ namespace MeshWorker Functional::assemble(const DoFInfo& info) { for (unsigned int i=0;i::assemble(const DoFInfo& info) { - for (unsigned int i=0;iblock(i)(info.face->user_index()) += info.J[i]; + results(1)->block(i)(info.face->user_index()) += info.value(i); else - results(0)->block(i)(info.cell->user_index()) += info.J[i]; + results(0)->block(i)(info.cell->user_index()) += info.value(i); } } @@ -1070,19 +1070,19 @@ namespace MeshWorker CellsAndFaces::assemble(const DoFInfo& info1, const DoFInfo& info2) { - for (unsigned int i=0;iblock(i)(info1.face->user_index()) += J; if (info2.face != info1.face) results(1)->block(i)(info2.face->user_index()) += J; } else { - results(0)->block(i)(info1.cell->user_index()) += .5*info1.J[i]; - results(0)->block(i)(info2.cell->user_index()) += .5*info2.J[i]; + results(0)->block(i)(info1.cell->user_index()) += .5*info1.value(i); + results(0)->block(i)(info2.cell->user_index()) += .5*info2.value(i); } } } @@ -1114,8 +1114,8 @@ namespace MeshWorker ResidualSimple::assemble(const DoFInfo& info) { for (unsigned int k=0;k& info) { for (unsigned int i=0;i::assemble(const DoFInfo& info) { - assemble(info.M1[0].matrix, info.indices, info.indices); + assemble(info.matrix(0,false).matrix, info.indices, info.indices); } @@ -1266,10 +1266,10 @@ namespace MeshWorker MatrixSimple::assemble(const DoFInfo& info1, const DoFInfo& info2) { - assemble(info1.M1[0].matrix, info1.indices, info1.indices); - assemble(info1.M2[0].matrix, info1.indices, info2.indices); - assemble(info2.M1[0].matrix, info2.indices, info2.indices); - assemble(info2.M2[0].matrix, info2.indices, info1.indices); + assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices); + assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices); + assemble(info2.matrix(0,false).matrix, info2.indices, info2.indices); + assemble(info2.matrix(0,true).matrix, info2.indices, info1.indices); } @@ -1354,7 +1354,7 @@ namespace MeshWorker MGMatrixSimple::assemble(const DoFInfo& info) { const unsigned int level = info.cell->level(); - assemble((*matrix)[level], info.M1[0].matrix, info.indices, info.indices); + assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices); } @@ -1369,10 +1369,10 @@ namespace MeshWorker if (level1 == level2) { - assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices); - assemble((*matrix)[level1], info1.M2[0].matrix, info1.indices, info2.indices); - assemble((*matrix)[level1], info2.M1[0].matrix, info2.indices, info2.indices); - assemble((*matrix)[level1], info2.M2[0].matrix, info2.indices, info1.indices); + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); + assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices); + assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices); + assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); } else { @@ -1380,9 +1380,9 @@ namespace MeshWorker // Do not add info2.M1, // which is done by // the coarser cell - assemble((*matrix)[level1], info1.M1[0].matrix, info1.indices, info1.indices); - assemble_transpose((*flux_up)[level1],info1.M2[0].matrix, info2.indices, info1.indices); - assemble((*flux_down)[level1], info2.M2[0].matrix, info2.indices, info1.indices); + assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices); + assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices); + assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices); } } @@ -1474,7 +1474,7 @@ namespace MeshWorker const unsigned int row = matrices[i]->row; const unsigned int col = matrices[i]->column; - assemble(matrices[i]->matrix, info.M1[i].matrix, row, col, info.indices, info.indices); + assemble(matrices[i]->matrix, info.matrix(i,false).matrix, row, col, info.indices, info.indices); } } @@ -1493,10 +1493,10 @@ namespace MeshWorker const unsigned int row = matrices[i]->row; const unsigned int col = matrices[i]->column; - assemble(matrices[i]->matrix, info1.M1[i].matrix, row, col, info1.indices, info1.indices); - assemble(matrices[i]->matrix, info1.M2[i].matrix, row, col, info1.indices, info2.indices); - assemble(matrices[i]->matrix, info2.M1[i].matrix, row, col, info2.indices, info2.indices); - assemble(matrices[i]->matrix, info2.M2[i].matrix, row, col, info2.indices, info1.indices); + assemble(matrices[i]->matrix, info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices); + assemble(matrices[i]->matrix, info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices); + assemble(matrices[i]->matrix, info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices); + assemble(matrices[i]->matrix, info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices); } } @@ -1607,7 +1607,7 @@ namespace MeshWorker const unsigned int row = matrices[i]->row; const unsigned int col = matrices[i]->column; - assemble(matrices[i]->matrix[level], info.M1[i].matrix, row, col, + assemble(matrices[i]->matrix[level], info.matrix(i,false).matrix, row, col, info.indices, info.indices, level, level); } } @@ -1632,10 +1632,10 @@ namespace MeshWorker if (level1 == level2) { - assemble(matrices[i]->matrix[level1], info1.M1[i].matrix, row, col, info1.indices, info1.indices, level1, level1); - assemble(matrices[i]->matrix[level1], info1.M2[i].matrix, row, col, info1.indices, info2.indices, level1, level2); - assemble(matrices[i]->matrix[level1], info2.M1[i].matrix, row, col, info2.indices, info2.indices, level2, level2); - assemble(matrices[i]->matrix[level1], info2.M2[i].matrix, row, col, info2.indices, info1.indices, level2, level1); + assemble(matrices[i]->matrix[level1], info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices, level1, level1); + assemble(matrices[i]->matrix[level1], info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices, level1, level2); + assemble(matrices[i]->matrix[level1], info2.matrix(i,false).matrix, row, col, info2.indices, info2.indices, level2, level2); + assemble(matrices[i]->matrix[level1], info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices, level2, level1); } else { @@ -1645,11 +1645,11 @@ namespace MeshWorker // Do not add M22, // which is done by // the coarser cell - assemble(matrices[i]->matrix[level1], info1.M1[i].matrix, row, col, + assemble(matrices[i]->matrix[level1], info1.matrix(i,false).matrix, row, col, info1.indices, info1.indices, level1, level1); - assemble(flux_up[i]->matrix[level1], info1.M2[i].matrix, row, col, + assemble(flux_up[i]->matrix[level1], info1.matrix(i,true).matrix, row, col, info1.indices, info2.indices, level1, level2, true); - assemble(flux_down[i]->matrix[level1], info2.M2[i].matrix, row, col, + assemble(flux_down[i]->matrix[level1], info2.matrix(i,true).matrix, row, col, info2.indices, info1.indices, level2, level1); } } 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 8028b762f4..96219c8ae7 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_info.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.h @@ -96,6 +96,66 @@ namespace MeshWorker */ void reinit(const BlockIndices& local_sizes); + /** + * The number of scalar values. + */ + unsigned int n_values () const; + + /** + * The number of vectors. + */ + unsigned int n_vectors () const; + + /** + * The number of matrices. + */ + unsigned int n_matrices () const; + + /** + * Access scalar value at index + * @p i. + */ + number& value(unsigned int i); + + /** + * Read scalar value at index + * @p i. + */ + number value(unsigned int i) const; + + /** + * Access vector at index @p i. + */ + BlockVector& vector(unsigned int i); + + /** + * Read vector at index @p i. + */ + const BlockVector& vector(unsigned int i) const; + + /** + * Access matrix at index @p + * i. For results on internal + * faces, a true value for @p + * external refers to the flux + * between cells, while false + * refers to entries coupling + * inside the cell. + */ + MatrixBlock >& matrix(unsigned int i, bool external = false); + + /** + * Read matrix at index @p + * i. For results on internal + * faces, a true value for @p + * external refers to the flux + * between cells, while false + * refers to entries coupling + * inside the cell. + */ + const MatrixBlock >& matrix(unsigned int i, bool external = false) const; + + private: /** * The local numbers, * computed on a cell or on a @@ -135,8 +195,8 @@ namespace MeshWorker /** - * Basic info class only containing information on geometry and - * degrees of freedom of the mesh object. + * A class containing information on geometry and degrees of freedom + * of a mesh object. * * The information in these objects is usually used by one of the * Assembler classes. It is also the kind of information which is @@ -145,15 +205,9 @@ namespace MeshWorker * * In addition to the information on degrees of freedom stored in this * class, it also provides the local computation space for the worker - * object operating on it. This space is provided by the base class - * template DATATYPE. This base class will automatically + * object operating on it in LocalResults. This base class will automatically * reinitialized on each cell, but initial setup is up to the user and - * should be done when initialize() for this class is called. The - * currently available base classes are - *
    - *
  • LocalVectors - *
  • LocalMatrices - *
+ * should be done when initialize() for this class is called. * * This class operates in two different modes, corresponding to the * data models discussed in the Assembler namespace documentation. @@ -163,7 +217,9 @@ namespace MeshWorker * BlockInfo::initialize_local(). If this function has been used, or * the vector has been changed from zero-length, then local dof * indices stored in this object will automatically be renumbered to - * reflect local block structure. + * reflect local block structure. This means, the first entries in + * #indices will refer to the first block of the system, then comes + * the second block and so on. * * The BlockInfo object is stored as a pointer. Therefore, if the * block structure changes, for instance because of mesh refinement, @@ -258,7 +314,7 @@ namespace MeshWorker SmartPointer > block_info; private: - /// Fill index vector + /// Fill index vector with active indices void get_indices(const typename DoFHandler::cell_iterator& c); /// Fill index vector with level indices @@ -329,24 +385,16 @@ namespace MeshWorker * @author Guido Kanschat, 2009 */ template - class IntegrationInfo : public DoFInfo + class IntegrationInfo { private: /// vector of FEValues objects std::vector > fevalv; public: /** - * Constructor forwarding - * information to DoFInfo. + * Constructor. */ - IntegrationInfo(const BlockInfo& block_info); - - /** - * Constructor forwarding - * information to DoFInfo. - */ - template - IntegrationInfo(const DH& dof_handler); + IntegrationInfo(); /** * Build all internal @@ -376,7 +424,8 @@ namespace MeshWorker void initialize(const FiniteElement& el, const Mapping& mapping, const Quadrature& quadrature, - const UpdateFlags flags); + const UpdateFlags flags, + const BlockInfo* local_block_info = 0); /** * Initialize the data @@ -467,41 +516,16 @@ namespace MeshWorker * Reinitialize internal data * structures for use on a cell. */ - template - void reinit(const DHCellIterator& c); - + void reinit(const DoFInfo& i); + /** - * Reinitialize internal data - * structures for use on a face. - */ - template - void reinit(const DHCellIterator& c, - const DHFaceIterator& f, - const unsigned int fn); - - /** - * Reinitialize internal data - * structures for use on a subface. - */ - template - void reinit(const DHCellIterator& c, - const DHFaceIterator& f, - const unsigned int fn, - const unsigned int sn); - - - /** - * @deprecated This is the - * old version not using - * VectorSelector. - * * Use the finite element * functions in #global_data * and fill the vectors * #values, #gradients and * #hessians. */ - void fill_local_data(const bool split_fevalues); + void fill_local_data(const DoFInfo& info, const bool split_fevalues); /** * The global data vector @@ -554,15 +578,13 @@ namespace MeshWorker template IntegrationInfoBox(const T&); - template + template void initialize(const WORKER&, - ASSEMBLER &assembler, const FiniteElement& el, const Mapping& mapping); - template + template void initialize(const WORKER&, - ASSEMBLER &assembler, const FiniteElement& el, const Mapping& mapping, const NamedData& data); @@ -573,6 +595,7 @@ namespace MeshWorker boost::shared_ptr > boundary_data; boost::shared_ptr > face_data; + DoFInfo dof_info; CellInfo cell_info; FaceInfo boundary_info; FaceInfo face_info; @@ -663,6 +686,101 @@ namespace MeshWorker } + template + inline + unsigned int + LocalResults::n_values() const + { + return J.size(); + } + + + template + inline + unsigned int + LocalResults::n_vectors() const + { + return R.size(); + } + + + template + inline + unsigned int + LocalResults::n_matrices() const + { + return M1.size(); + } + + + template + inline + number& + LocalResults::value(unsigned int i) + { + AssertIndexRange(i,J.size()); + return J[i]; + } + + + template + inline + BlockVector& + LocalResults::vector(unsigned int i) + { + AssertIndexRange(i,R.size()); + return R[i]; + } + + + template + inline + MatrixBlock >& + LocalResults::matrix(unsigned int i, bool external) + { + if (external) + { + AssertIndexRange(i,M2.size()); + return M2[i]; + } + AssertIndexRange(i,M1.size()); + return M1[i]; + } + + template + inline + number + LocalResults::value(unsigned int i) const + { + AssertIndexRange(i,J.size()); + return J[i]; + } + + + template + inline + const BlockVector& + LocalResults::vector(unsigned int i) const + { + AssertIndexRange(i,R.size()); + return R[i]; + } + + + template + inline + const MatrixBlock >& + LocalResults::matrix(unsigned int i, bool external) const + { + if (external) + { + AssertIndexRange(i,M2.size()); + return M2[i]; + } + AssertIndexRange(i,M1.size()); + return M1[i]; + } + //----------------------------------------------------------------------// template @@ -746,17 +864,6 @@ namespace MeshWorker //----------------------------------------------------------------------// - template - template - IntegrationInfo::IntegrationInfo(const DH& dof_handler) - : - DoFInfo(dof_handler), - fevalv(0), - multigrid(false), - global_data(boost::shared_ptr >(new VectorDataBase)) - {} - - template inline const FVB& IntegrationInfo::fe_values() const @@ -776,94 +883,55 @@ namespace MeshWorker template - template - inline void - IntegrationInfo::reinit(const DHCellIterator& c) - { - DoFInfo::reinit(c); - for (unsigned int i=0;i& fe = dynamic_cast&> (febase); - fe.reinit(this->cell); - } - - const bool split_fevalues = this->block_info != 0; - fill_local_data(split_fevalues); - } - - - template - template inline void - IntegrationInfo::reinit( - const DHCellIterator& c, - const DHFaceIterator& f, - const unsigned int fn) + IntegrationInfo::reinit(const DoFInfo& info) { - DoFInfo::reinit(c, f, fn); for (unsigned int i=0;i& fe = dynamic_cast&> (febase); - fe.reinit(this->cell, fn); + if (info.sub_number != deal_II_numbers::invalid_unsigned_int) + { + // This is a subface + FESubfaceValues& fe = dynamic_cast&> (febase); + fe.reinit(info.cell, info.face_number, info.sub_number); + } + else if (info.face_number != deal_II_numbers::invalid_unsigned_int) + { + // This is a face + FEFaceValues& fe = dynamic_cast&> (febase); + fe.reinit(info.cell, info.face_number); + } + else + { + // This is a cell + FEValues& fe = dynamic_cast&> (febase); + fe.reinit(info.cell); + } } - const bool split_fevalues = this->block_info != 0; - fill_local_data(split_fevalues); + const bool split_fevalues = info.block_info != 0; + fill_local_data(info, split_fevalues); } - template - template - inline void - IntegrationInfo::reinit( - const DHCellIterator& c, - const DHFaceIterator& f, - const unsigned int fn, - const unsigned int sn) - { - DoFInfo::reinit(c, f, fn, sn); - for (unsigned int i=0;i& fe = dynamic_cast&> (febase); - fe.reinit(this->cell, fn, sn); - } - - const bool split_fevalues = this->block_info != 0; - fill_local_data(split_fevalues); - } - //----------------------------------------------------------------------// template template IntegrationInfoBox::IntegrationInfoBox(const T& t) : - cell_info(t), - boundary_info(t), - face_info(t), - subface_info(t), - neighbor_info(t) + dof_info(t) {} template - template + template void IntegrationInfoBox:: initialize(const WORKER& integrator, - ASSEMBLER &assembler, const FiniteElement& el, const Mapping& mapping) { - assembler.initialize_info(cell_info, false); - assembler.initialize_info(boundary_info, false); - assembler.initialize_info(face_info, true); - assembler.initialize_info(subface_info, true); - assembler.initialize_info(neighbor_info, true); - cell_info.initialize >(el, mapping, integrator.cell_quadrature, integrator.cell_flags); boundary_info.initialize >(el, mapping, integrator.boundary_quadrature, @@ -878,16 +946,15 @@ namespace MeshWorker template - template + template void IntegrationInfoBox::initialize( const WORKER& integrator, - ASSEMBLER &assembler, const FiniteElement& el, const Mapping& mapping, const NamedData& data) { - initialize(integrator, assembler, el, mapping); + initialize(integrator, el, mapping); boost::shared_ptr > p; p = boost::shared_ptr >(new VectorData (integrator.cell_selector)); diff --git a/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h b/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h index 3061cff461..82baecbaa8 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_info.templates.h @@ -63,9 +63,8 @@ namespace MeshWorker //----------------------------------------------------------------------// template - IntegrationInfo::IntegrationInfo(const BlockInfo& block_info) + IntegrationInfo::IntegrationInfo() : - DoFInfo(block_info), fevalv(0), multigrid(false), global_data(boost::shared_ptr >(new VectorDataBase)) @@ -79,9 +78,10 @@ namespace MeshWorker const FiniteElement& el, const Mapping& mapping, const Quadrature& quadrature, - const UpdateFlags flags) + const UpdateFlags flags, + const BlockInfo* block_info) { - if (this->block_info == 0 || this->block_info->local().size() == 0) + if (block_info == 0 || block_info->local().size() == 0) { fevalv.resize(1); fevalv[0] = boost::shared_ptr ( @@ -165,21 +165,21 @@ namespace MeshWorker template void - IntegrationInfo::fill_local_data(bool split_fevalues) + IntegrationInfo::fill_local_data(const DoFInfo& info, bool split_fevalues) { if (split_fevalues) { unsigned int comp = 0; // Loop over all blocks - for (unsigned int b=0;bblock_info->local().size();++b) + for (unsigned int b=0;blocal().size();++b) { - const unsigned int fe_no = this->block_info->base_element(b); + const unsigned int fe_no = info.block_info->base_element(b); const FEValuesBase& fe = this->fe_values(fe_no); const unsigned int n_comp = fe.get_fe().n_components(); - const unsigned int block_start = this->block_info->local().block_start(b); - const unsigned int block_size = this->block_info->local().block_size(b); + const unsigned int block_start = info.block_info->local().block_start(b); + const unsigned int block_size = info.block_info->local().block_size(b); - this->global_data->fill(values, gradients, hessians, fe, this->indices, + this->global_data->fill(values, gradients, hessians, fe, info.indices, comp, n_comp, block_start, block_size); comp += n_comp; } @@ -188,8 +188,8 @@ namespace MeshWorker { const FEValuesBase& fe = this->fe_values(0); const unsigned int n_comp = fe.get_fe().n_components(); - this->global_data->fill(values, gradients, hessians, fe, this->indices, - 0, n_comp, 0, this->indices.size()); + this->global_data->fill(values, gradients, hessians, fe, info.indices, + 0, n_comp, 0, info.indices.size()); } } } diff --git a/deal.II/deal.II/include/numerics/mesh_worker_loop.h b/deal.II/deal.II/include/numerics/mesh_worker_loop.h index 140ac08a11..b48f92cc37 100644 --- a/deal.II/deal.II/include/numerics/mesh_worker_loop.h +++ b/deal.II/deal.II/include/numerics/mesh_worker_loop.h @@ -66,14 +66,15 @@ namespace MeshWorker * @ingroup MeshWorker * @author Guido Kanschat, 2009 */ - template + template void loop(ITERATOR begin, typename identity::type end, + DOFINFO cell_dof_info, CELLINFO& cell_info, FACEINFO& boundary_info, FACEINFO& face_info, FACEINFO& subface_info, FACEINFO& neighbor_info, - const std_cxx1x::function &cell_worker, - const std_cxx1x::function &boundary_worker, - const std_cxx1x::function &face_worker, + const std_cxx1x::function &cell_worker, + const std_cxx1x::function &boundary_worker, + const std_cxx1x::function &face_worker, ASSEMBLER &assembler, bool cells_first = true) { @@ -81,6 +82,13 @@ namespace MeshWorker const bool integrate_boundary = (boundary_worker != 0); const bool integrate_interior_face = (face_worker != 0); + ; + DOFINFO face_dof_info = cell_dof_info; + DOFINFO neighbor_dof_info = cell_dof_info; + assembler.initialize_info(cell_dof_info, false); + assembler.initialize_info(face_dof_info, true); + assembler.initialize_info(neighbor_dof_info, true); + // Loop over all cells for (ITERATOR cell = begin; cell != end; ++cell) { @@ -89,9 +97,10 @@ namespace MeshWorker // before faces if (integrate_cell && cells_first) { - cell_info.reinit(cell); - cell_worker(cell_info); - assembler.assemble (cell_info); + cell_dof_info.reinit(cell); + cell_info.reinit(cell_dof_info); + cell_worker(cell_dof_info, cell_info); + assembler.assemble(cell_dof_info); } if (integrate_interior_face || integrate_boundary) @@ -102,9 +111,10 @@ namespace MeshWorker { if (integrate_boundary) { - boundary_info.reinit(cell, face, face_no); - boundary_worker(boundary_info); - assembler.assemble (boundary_info); + cell_dof_info.reinit(cell, face, face_no); + boundary_info.reinit(cell_dof_info); + boundary_worker(cell_dof_info, boundary_info); + assembler.assemble(cell_dof_info); } } else if (integrate_interior_face) @@ -132,14 +142,18 @@ namespace MeshWorker = cell->neighbor_of_coarser_neighbor(face_no); typename ITERATOR::AccessorType::Container::face_iterator nface = neighbor->face(neighbor_face_no.first); - face_info.reinit(cell, face, face_no); - subface_info.reinit(neighbor, nface, neighbor_face_no.first, neighbor_face_no.second); + + face_dof_info.reinit(cell, face, face_no); + face_info.reinit(face_dof_info); + neighbor_dof_info.reinit(neighbor, nface, + neighbor_face_no.first, neighbor_face_no.second); + subface_info.reinit(neighbor_dof_info); // Neighbor // first to // conform to // old version - face_worker(face_info, subface_info); - assembler.assemble (face_info, subface_info); + face_worker(face_dof_info, neighbor_dof_info, face_info, subface_info); + assembler.assemble (face_dof_info, neighbor_dof_info); } else { @@ -162,12 +176,13 @@ namespace MeshWorker unsigned int neighbor_face_no = cell->neighbor_of_neighbor(face_no); Assert (neighbor->face(neighbor_face_no) == face, ExcInternalError()); // Regular interior face - face_info.reinit(cell, face, face_no); - neighbor_info.reinit(neighbor, neighbor->face(neighbor_face_no), - neighbor_face_no); - - face_worker(face_info, neighbor_info); - assembler.assemble (face_info, neighbor_info); + face_dof_info.reinit(cell, face, face_no); + face_info.reinit(face_dof_info); + neighbor_dof_info.reinit(neighbor, neighbor->face(neighbor_face_no), + neighbor_face_no); + neighbor_info.reinit(neighbor_dof_info); + face_worker(face_dof_info, neighbor_dof_info, face_info, neighbor_info); + assembler.assemble(face_dof_info, neighbor_dof_info); } } } // faces @@ -175,9 +190,10 @@ namespace MeshWorker // have to be handled first if (integrate_cell && !cells_first) { - cell_info.reinit(cell); - cell_worker(cell_info); - assembler.assemble (cell_info); + cell_dof_info.reinit(cell); + cell_info.reinit(cell_dof_info); + cell_worker(cell_dof_info, cell_info); + assembler.assemble (cell_dof_info); } } } @@ -192,14 +208,15 @@ namespace MeshWorker void integration_loop(ITERATOR begin, typename identity::type end, IntegrationInfoBox& box, - const std_cxx1x::function &cell_worker, - const std_cxx1x::function &boundary_worker, - const std_cxx1x::function &face_worker, + const std_cxx1x::function&, CELLINFO &)> &cell_worker, + const std_cxx1x::function&, FACEINFO &)> &boundary_worker, + const std_cxx1x::function&, DoFInfo&, FACEINFO &, FACEINFO &)> &face_worker, ASSEMBLER &assembler, bool cells_first = true) { - loop + loop,CELLINFO,FACEINFO> (begin, end, + box.dof_info, box.cell_info, box.boundary_info, box.face_info, box.subface_info, box.neighbor_info, diff --git a/deal.II/deal.II/source/numerics/mesh_worker_info.cc b/deal.II/deal.II/source/numerics/mesh_worker_info.cc index cdd5bb741c..9d418fd23b 100644 --- a/deal.II/deal.II/source/numerics/mesh_worker_info.cc +++ b/deal.II/deal.II/source/numerics/mesh_worker_info.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 2009 by the deal.II authors +// Copyright (C) 2009, 2010 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -34,15 +34,15 @@ namespace MeshWorker template void IntegrationInfo > ::initialize >( const FiniteElement&, const Mapping&, - const Quadrature::integral_dimension>&, const UpdateFlags); + const Quadrature::integral_dimension>&, const UpdateFlags, const BlockInfo*); template void IntegrationInfo > ::initialize >( const FiniteElement&, const Mapping&, - const Quadrature::integral_dimension>&, const UpdateFlags); + const Quadrature::integral_dimension>&, const UpdateFlags, const BlockInfo*); template void IntegrationInfo > ::initialize >( const FiniteElement&, const Mapping&, - const Quadrature::integral_dimension>&, const UpdateFlags); + const Quadrature::integral_dimension>&, const UpdateFlags, const BlockInfo*); } #endif diff --git a/deal.II/examples/step-38/step-38.cc b/deal.II/examples/step-38/step-38.cc index ab3f260722..9d7a1ad814 100644 --- a/deal.II/examples/step-38/step-38.cc +++ b/deal.II/examples/step-38/step-38.cc @@ -39,6 +39,7 @@ // Here come the new include files // for using the MeshWorker framework: #include +#include #include #include @@ -200,9 +201,10 @@ class DGMethod // types of arguments, but have // in fact other arguments // already bound. - static void integrate_cell_term (CellInfo& info); - static void integrate_boundary_term (FaceInfo& info); - static void integrate_face_term (FaceInfo& info1, + static void integrate_cell_term (MeshWorker::DoFInfo& dinfo, CellInfo& info); + static void integrate_boundary_term (MeshWorker::DoFInfo& dinfo, FaceInfo& info); + static void integrate_face_term (MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, FaceInfo& info1, FaceInfo& info2); }; @@ -378,19 +380,20 @@ void DGMethod::assemble_system () // receives all the stuff we // created so far. MeshWorker::IntegrationInfoBox info_box(dof_handler); - info_box.initialize(integration_worker, assembler, fe, mapping); + info_box.initialize(integration_worker, fe, mapping); // Finally, the integration loop // over all active cells // (determined by the first - // argument, which is an active iterator). - MeshWorker::integration_loop + // argument, which is an active + // iterator). + MeshWorker::integration_loop (dof_handler.begin_active(), dof_handler.end(), info_box, &DGMethod::integrate_cell_term, &DGMethod::integrate_boundary_term, &DGMethod::integrate_face_term, - assembler); + assembler, true); } @@ -412,7 +415,7 @@ void DGMethod::assemble_system () // added soon). template -void DGMethod::integrate_cell_term (CellInfo& info) +void DGMethod::integrate_cell_term (MeshWorker::DoFInfo& dinfo, CellInfo& info) { // First, let us retrieve some of // the objects used here from @@ -422,7 +425,7 @@ void DGMethod::integrate_cell_term (CellInfo& info) // looks more complicated than // might seem necessary. const FEValuesBase& fe_v = info.fe_values(); - FullMatrix& local_matrix = info.M1[0].matrix; + FullMatrix& local_matrix = dinfo.matrix(0).matrix; const std::vector &JxW = fe_v.get_JxW_values (); // With these objects, we continue @@ -458,11 +461,11 @@ void DGMethod::integrate_cell_term (CellInfo& info) // FESubfaceValues, in order to get access to // normal vectors. template -void DGMethod::integrate_boundary_term (FaceInfo& info) +void DGMethod::integrate_boundary_term (MeshWorker::DoFInfo& dinfo, FaceInfo& info) { const FEFaceValuesBase& fe_v = info.fe_values(); - FullMatrix& local_matrix = info.M1[0].matrix; - Vector& local_vector = info.R[0].block(0); + FullMatrix& local_matrix = dinfo.matrix(0).matrix; + Vector& local_vector = dinfo.vector(0).block(0); const std::vector &JxW = fe_v.get_JxW_values (); const std::vector > &normals = fe_v.get_normal_vectors (); @@ -504,7 +507,9 @@ void DGMethod::integrate_boundary_term (FaceInfo& info) // for each cell and two for coupling // back and forth. template -void DGMethod::integrate_face_term (FaceInfo& info1, +void DGMethod::integrate_face_term (MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + FaceInfo& info1, FaceInfo& info2) { // For quadrature points, weights, @@ -531,10 +536,10 @@ void DGMethod::integrate_face_term (FaceInfo& info1, // interior couplings of that cell, // while the second contains the // couplings between cells. - FullMatrix& u1_v1_matrix = info1.M1[0].matrix; - FullMatrix& u2_v1_matrix = info1.M2[0].matrix; - FullMatrix& u1_v2_matrix = info2.M2[0].matrix; - FullMatrix& u2_v2_matrix = info2.M1[0].matrix; + FullMatrix& u1_v1_matrix = dinfo1.matrix(0,false).matrix; + FullMatrix& u2_v1_matrix = dinfo1.matrix(0,true).matrix; + FullMatrix& u1_v2_matrix = dinfo2.matrix(0,true).matrix; + FullMatrix& u2_v2_matrix = dinfo2.matrix(0,false).matrix; // Here, following the previous // functions, we would have the diff --git a/deal.II/examples/step-39/step-39.cc b/deal.II/examples/step-39/step-39.cc index c39e0a8b74..7de4081c58 100644 --- a/deal.II/examples/step-39/step-39.cc +++ b/deal.II/examples/step-39/step-39.cc @@ -97,9 +97,13 @@ template class MatrixIntegrator : public Subscriptor { public: - static void cell(typename MeshWorker::IntegrationWorker::CellInfo& info); - static void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info); - static void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + static void cell(MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationWorker::CellInfo& info); + static void bdry(MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationWorker::FaceInfo& info); + static void face(MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + typename MeshWorker::IntegrationWorker::FaceInfo& info1, typename MeshWorker::IntegrationWorker::FaceInfo& info2); }; @@ -112,10 +116,12 @@ class MatrixIntegrator : public Subscriptor // degrees of freedom associated // with the shape functions. template -void MatrixIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) +void MatrixIntegrator::cell( + MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationWorker::CellInfo& info) { const FEValuesBase& fe = info.fe_values(); - FullMatrix& local_matrix = info.M1[0].matrix; + FullMatrix& local_matrix = dinfo.matrix(0,false).matrix; for (unsigned int k=0; k::cell(typename MeshWorker::IntegrationWorker::Ce // On boundary faces, we use the // Nitsche boundary condition template -void MatrixIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) +void MatrixIntegrator::bdry( + MeshWorker::DoFInfo& dinfo, + typename MeshWorker::IntegrationWorker::FaceInfo& info) { const FEFaceValuesBase& fe = info.fe_values(); - FullMatrix& local_matrix = info.M1[0].matrix; + FullMatrix& local_matrix = dinfo.matrix(0,false).matrix; const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); for (unsigned k=0;k::bdry(typename MeshWorker::IntegrationWorker::Fa template -void MatrixIntegrator::face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, - typename MeshWorker::IntegrationWorker::FaceInfo& info2) +void MatrixIntegrator::face( + MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + typename MeshWorker::IntegrationWorker::FaceInfo& info1, + typename MeshWorker::IntegrationWorker::FaceInfo& info2) { const FEFaceValuesBase& fe1 = info1.fe_values(); const FEFaceValuesBase& fe2 = info2.fe_values(); - FullMatrix& matrix_v1u1 = info1.M1[0].matrix; - FullMatrix& matrix_v1u2 = info1.M2[0].matrix; - FullMatrix& matrix_v2u1 = info2.M2[0].matrix; - FullMatrix& matrix_v2u2 = info2.M1[0].matrix; + FullMatrix& matrix_v1u1 = dinfo1.matrix(0,false).matrix; + FullMatrix& matrix_v1u2 = dinfo1.matrix(0,true).matrix; + FullMatrix& matrix_v2u1 = dinfo2.matrix(0,true).matrix; + FullMatrix& matrix_v2u2 = dinfo2.matrix(0,false).matrix; const unsigned int deg = fe1.get_fe().tensor_degree(); - double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure(); - double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure(); - if (info1.cell->has_children() ^ info2.cell->has_children()) + double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); + if (dinfo1.cell->has_children() ^ dinfo2.cell->has_children()) { - Assert (info1.face == info2.face, ExcInternalError()); - Assert (info1.face->has_children(), ExcInternalError()); -// Assert (info1.cell->has_children(), ExcInternalError()); + Assert (dinfo1.face == dinfo2.face, ExcInternalError()); + Assert (dinfo1.face->has_children(), ExcInternalError()); +// Assert (dinfo1.cell->has_children(), ExcInternalError()); penalty1 *= 2; } const double penalty = penalty1 + penalty2; @@ -196,29 +207,31 @@ template class RHSIntegrator : public Subscriptor { public: - static void cell(typename MeshWorker::IntegrationWorker::CellInfo& info); - static void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info); - static void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + static void cell(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationWorker::CellInfo& info); + static void bdry(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationWorker::FaceInfo& info); + static void face(MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + typename MeshWorker::IntegrationWorker::FaceInfo& info1, typename MeshWorker::IntegrationWorker::FaceInfo& info2); }; template -void RHSIntegrator::cell(typename MeshWorker::IntegrationWorker::CellInfo&) +void RHSIntegrator::cell(MeshWorker::DoFInfo&, typename MeshWorker::IntegrationWorker::CellInfo&) {} template -void RHSIntegrator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) +void RHSIntegrator::bdry(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationWorker::FaceInfo& info) { const FEFaceValuesBase& fe = info.fe_values(); - Vector& local_vector = info.R[0].block(0); + Vector& local_vector = dinfo.vector(0).block(0); std::vector boundary_values(fe.n_quadrature_points); exact_solution.value_list(fe.get_quadrature_points(), boundary_values); const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); for (unsigned k=0;k::bdry(typename MeshWorker::IntegrationWorker::FaceI template -void RHSIntegrator::face(typename MeshWorker::IntegrationWorker::FaceInfo&, +void RHSIntegrator::face(MeshWorker::DoFInfo&, + MeshWorker::DoFInfo&, + typename MeshWorker::IntegrationWorker::FaceInfo&, typename MeshWorker::IntegrationWorker::FaceInfo&) {} @@ -238,29 +253,31 @@ template class Estimator : public Subscriptor { public: - static void cell(typename MeshWorker::IntegrationWorker::CellInfo& info); - static void bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info); - static void face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, + static void cell(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationWorker::CellInfo& info); + static void bdry(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationWorker::FaceInfo& info); + static void face(MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + typename MeshWorker::IntegrationWorker::FaceInfo& info1, typename MeshWorker::IntegrationWorker::FaceInfo& info2); }; template -void Estimator::cell(typename MeshWorker::IntegrationWorker::CellInfo& info) +void Estimator::cell(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationWorker::CellInfo& info) { const FEValuesBase& fe = info.fe_values(); const std::vector >& DDuh = info.hessians[0][0]; for (unsigned k=0;kdiameter() * trace(DDuh[k]); - info.J[0] += t*t * fe.JxW(k); + const double t = dinfo.cell->diameter() * trace(DDuh[k]); + dinfo.value(0) += t*t * fe.JxW(k); } } template -void Estimator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& info) +void Estimator::bdry(MeshWorker::DoFInfo& dinfo, typename MeshWorker::IntegrationWorker::FaceInfo& info) { const FEFaceValuesBase& fe = info.fe_values(); @@ -270,16 +287,18 @@ void Estimator::bdry(typename MeshWorker::IntegrationWorker::FaceInfo& const std::vector& uh = info.values[0][0]; const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty = 2. * deg * (deg+1) * info.face->measure() / info.cell->measure(); + const double penalty = 2. * deg * (deg+1) * dinfo.face->measure() / dinfo.cell->measure(); for (unsigned k=0;k -void Estimator::face(typename MeshWorker::IntegrationWorker::FaceInfo& info1, +void Estimator::face(MeshWorker::DoFInfo& dinfo1, + MeshWorker::DoFInfo& dinfo2, + typename MeshWorker::IntegrationWorker::FaceInfo& info1, typename MeshWorker::IntegrationWorker::FaceInfo& info2) { const FEFaceValuesBase& fe = info1.fe_values(); @@ -289,19 +308,19 @@ void Estimator::face(typename MeshWorker::IntegrationWorker::FaceInfo& const std::vector >& Duh2 = info2.gradients[0][0]; const unsigned int deg = fe.get_fe().tensor_degree(); - const double penalty1 = deg * (deg+1) * info1.face->measure() / info1.cell->measure(); - const double penalty2 = deg * (deg+1) * info2.face->measure() / info2.cell->measure(); + const double penalty1 = deg * (deg+1) * dinfo1.face->measure() / dinfo1.cell->measure(); + const double penalty2 = deg * (deg+1) * dinfo2.face->measure() / dinfo2.cell->measure(); const double penalty = penalty1 + penalty2; - const double h = info1.face->measure(); + const double h = dinfo1.face->measure(); for (unsigned k=0;k::assemble_matrix() assembler.initialize(matrix); MeshWorker::IntegrationInfoBox info_box(dof_handler); - info_box.initialize(integration_worker, assembler, fe, mapping); - MeshWorker::integration_loop( + info_box.initialize(integration_worker, fe, mapping); + MeshWorker::integration_loop( dof_handler.begin_active(), dof_handler.end(), info_box, &MatrixIntegrator::cell, @@ -451,8 +470,8 @@ Step39::assemble_mg_matrix() assembler.initialize(mg_matrix); assembler.initialize_fluxes(mg_matrix_dg_up, mg_matrix_dg_down); MeshWorker::IntegrationInfoBox info_box(mg_dof_handler); - info_box.initialize(integration_worker, assembler, fe, mapping); - MeshWorker::integration_loop( + info_box.initialize(integration_worker, fe, mapping); + MeshWorker::integration_loop( mg_dof_handler.begin(), mg_dof_handler.end(), info_box, &MatrixIntegrator::cell, @@ -479,9 +498,9 @@ Step39::assemble_right_hand_side() data.add(rhs, "RHS"); assembler.initialize(data); MeshWorker::IntegrationInfoBox info_box(dof_handler); - info_box.initialize(integration_worker, assembler, fe, mapping); + info_box.initialize(integration_worker, fe, mapping); - MeshWorker::integration_loop( + MeshWorker::integration_loop( dof_handler.begin_active(), dof_handler.end(), info_box, &RHSIntegrator::cell, @@ -600,8 +619,8 @@ Step39::estimate() out_data.add(est, "cells"); assembler.initialize(out_data, false); MeshWorker::IntegrationInfoBox info_box(dof_handler); - info_box.initialize(integration_worker, assembler, fe, mapping, solution_data); - MeshWorker::integration_loop ( + info_box.initialize(integration_worker, fe, mapping, solution_data); + MeshWorker::integration_loop ( dof_handler.begin_active(), dof_handler.end(), info_box, &Estimator::cell, -- 2.39.5