From 1e603742a4c3c134c18f16e63f24052870c201f3 Mon Sep 17 00:00:00 2001 From: Baerbel Jannsen Date: Thu, 1 Jul 2010 12:40:12 +0000 Subject: [PATCH] include multigrid assembling for continuous elements to the MeshWorker git-svn-id: https://svn.dealii.org/trunk@21429 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/multigrid/mg_constraints.h | 82 +++- .../include/numerics/mesh_worker_assembler.h | 386 ++++++++++++------ 2 files changed, 326 insertions(+), 142 deletions(-) diff --git a/deal.II/deal.II/include/multigrid/mg_constraints.h b/deal.II/deal.II/include/multigrid/mg_constraints.h index c23d3e0ee0..35e46b12bc 100644 --- a/deal.II/deal.II/include/multigrid/mg_constraints.h +++ b/deal.II/deal.II/include/multigrid/mg_constraints.h @@ -16,6 +16,8 @@ #include #include +#include + #include #include @@ -68,7 +70,7 @@ class MGConstraints : public Subscriptor * is subject to a boundary * constraint. */ - bool at_boundary(unsigned int level, unsigned int index) const; + bool is_boundary_index(unsigned int level, unsigned int index) const; /** * Determine whether a dof index @@ -76,13 +78,22 @@ class MGConstraints : public Subscriptor */ bool at_refinement_edge(unsigned int level, unsigned int index) const; - private: + /** + * Determine whether a dof index + * is at the refinement edge and + * subject to a boundary + * constraint . + */ + bool at_refinement_edge_boundary(unsigned int level, unsigned int index) const; + /** * The indices of boundary dofs * for each level. */ std::vector > boundary_indices; + private: + /** * The degrees of freedom on the * refinement edge between a @@ -102,15 +113,74 @@ class MGConstraints : public Subscriptor std::vector > refinement_edge_boundary_indices; }; +template +void MGConstraints::initialize(const MGDoFHandler& dof) +{ + MGTools::extract_inner_interface_dofs (dof, refinement_edge_indices, + refinement_edge_boundary_indices); +} + +template +void MGConstraints::initialize(const MGDoFHandler& dof, + const typename FunctionMap::type& function_map, + const std::vector& component_mask) +{ + const unsigned int nlevels = dof.get_tria().n_levels(); + boundary_indices.resize(nlevels); + refinement_edge_indices.resize(nlevels); + refinement_edge_boundary_indices.resize(nlevels); + + for(unsigned int l=0; l #include #include -#include #include +#include +#include +#include DEAL_II_NAMESPACE_OPEN @@ -120,16 +121,16 @@ namespace MeshWorker * Assemble the local values * into the global vectors. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble both local values * into the global vectors. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); /** * The value of the ith entry @@ -217,24 +218,24 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, + template + void initialize_info(DOFINFO& info, bool interior_face) const; /** * Assemble the local values * into the global vectors. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble both local values * into the global vectors. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); /** * The value of the ith entry @@ -286,24 +287,24 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, + template + void initialize_info(DOFINFO& info, bool interior_face) const; /** * Assemble the local residuals * into the global residuals. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble both local residuals * into the global residuals. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); private: /** * The global residal vectors @@ -369,8 +370,8 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, + template + void initialize_info(DOFINFO& info, bool interior_face) const; @@ -378,16 +379,16 @@ namespace MeshWorker * Assemble the local residuals * into the global residuals. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble both local residuals * into the global residuals. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); private: /** * Assemble a single local @@ -472,8 +473,8 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, + template + void initialize_info(DOFINFO& info, bool interior_face) const; /** @@ -481,8 +482,8 @@ namespace MeshWorker * DoFInfo::M1[0] * into the global matrix. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble both local @@ -490,9 +491,9 @@ namespace MeshWorker * objects into the global * matrix. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); private: /** * Assemble a single matrix @@ -555,6 +556,12 @@ namespace MeshWorker */ void initialize(MGLevelObject& m); + /** + * Initialize the multilevel + * constraints. + */ + void initialize(const MGConstraints& mg_constraints); + /** * Initialize the matrices * #flux_up and #flux_down @@ -592,8 +599,8 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, + template + void initialize_info(DOFINFO& info, bool interior_face) const; /** @@ -601,8 +608,8 @@ namespace MeshWorker * DoFInfo::M1[0] * into the global matrix. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble both local @@ -610,9 +617,9 @@ namespace MeshWorker * objects into the global * matrices. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); private: /** * Assemble a single matrix @@ -621,7 +628,8 @@ namespace MeshWorker void assemble(MATRIX& G, const FullMatrix& M, const std::vector& i1, - const std::vector& i2); + const std::vector& i2, + const unsigned int level = -1); /** * Assemble a single matrix @@ -630,7 +638,18 @@ namespace MeshWorker void assemble_transpose(MATRIX& G, const FullMatrix& M, const std::vector& i1, - const std::vector& i2); + const std::vector& i2, + const unsigned int level = -1); + + /** + * Assemble a single matrix + * into a global matrix. + */ + void assemble_down(MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level = -1); /** * The global matrix being @@ -671,6 +690,10 @@ namespace MeshWorker * fine to coarse. */ SmartPointer,MGMatrixSimple > interface_down; + /** + * A pointer to the object containing constraints. + */ + SmartPointer > mg_constraints; /** * The smallest positive @@ -756,8 +779,8 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, + template + void initialize_info(DOFINFO& info, bool interior_face) const; @@ -765,16 +788,16 @@ namespace MeshWorker * Assemble the local matrices * into the global matrices. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble all local matrices * into the global matrices. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); private: /** @@ -899,8 +922,8 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, + template + void initialize_info(DOFINFO& info, bool interior_face) const; @@ -908,16 +931,16 @@ namespace MeshWorker * Assemble the local matrices * into the global matrices. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble all local matrices * into the global matrices. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); private: /** @@ -1034,16 +1057,16 @@ namespace MeshWorker * require additional data * being initialized. */ - template - void initialize_info(DoFInfo& info, bool) const; + template + void initialize_info(DOFINFO& info, bool) const; /** * Assemble the matrix * DoFInfo::M1[0] * into the global matrix. */ - template - void assemble(const DoFInfo& info); + template + void assemble(const DOFINFO& info); /** * Assemble both local @@ -1051,9 +1074,9 @@ namespace MeshWorker * objects into the global * matrix. */ - template - void assemble(const DoFInfo& info1, - const DoFInfo& info2); + template + void assemble(const DOFINFO& info1, + const DOFINFO& info2); }; @@ -1068,9 +1091,9 @@ namespace MeshWorker template - template + template inline void - Functional::assemble(const DoFInfo& info) + Functional::assemble(const DOFINFO& info) { for (unsigned int i=0;i - template + template inline void - Functional::assemble(const DoFInfo& info1, - const DoFInfo& info2) + Functional::assemble(const DOFINFO& info1, + const DOFINFO& info2) { for (unsigned int i=0;i - template + template inline void - CellsAndFaces::initialize_info(DoFInfo& info, bool) const + CellsAndFaces::initialize_info(DOFINFO& info, bool) const { info.initialize_numbers(results(0)->n_blocks()); } template - template + template inline void - CellsAndFaces::assemble(const DoFInfo& info) + CellsAndFaces::assemble(const DOFINFO& info) { for (unsigned int i=0;i - template + template inline void - CellsAndFaces::assemble(const DoFInfo& info1, - const DoFInfo& info2) + CellsAndFaces::assemble(const DOFINFO& info1, + const DOFINFO& info2) { for (unsigned int i=0;i - template + template inline void - ResidualSimple::initialize_info(DoFInfo& info, bool) const + ResidualSimple::initialize_info(DOFINFO& info, bool) const { info.initialize_vectors(residuals.size()); } template - template + template inline void - ResidualSimple::assemble(const DoFInfo& info) + ResidualSimple::assemble(const DOFINFO& info) { for (unsigned int k=0;k - template + template inline void - ResidualSimple::assemble(const DoFInfo& info1, - const DoFInfo& info2) + ResidualSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) { for (unsigned int k=0;k - template + template inline void ResidualLocalBlocksToGlobalBlocks::initialize_info( - DoFInfo& info, bool) const + DOFINFO& info, bool) const { info.initialize_vectors(residuals.size()); } @@ -1297,10 +1320,10 @@ namespace MeshWorker template - template + template inline void ResidualLocalBlocksToGlobalBlocks::assemble( - const DoFInfo& info) + const DOFINFO& info) { for (unsigned int i=0;i - template + template inline void ResidualLocalBlocksToGlobalBlocks::assemble( - const DoFInfo& info1, - const DoFInfo& info2) + const DOFINFO& info1, + const DOFINFO& info2) { for (unsigned int i=0;i - template + template inline void - MatrixSimple::initialize_info(DoFInfo& info, + MatrixSimple::initialize_info(DOFINFO& info, bool interior_face) const { info.initialize_matrices(1, interior_face); @@ -1382,19 +1405,19 @@ namespace MeshWorker template - template + template inline void - MatrixSimple::assemble(const DoFInfo& info) + MatrixSimple::assemble(const DOFINFO& info) { assemble(info.matrix(0,false).matrix, info.indices, info.indices); } template - template + template inline void - MatrixSimple::assemble(const DoFInfo& info1, - const DoFInfo& info2) + MatrixSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) { assemble(info1.matrix(0,false).matrix, info1.indices, info1.indices); assemble(info1.matrix(0,true).matrix, info1.indices, info2.indices); @@ -1420,6 +1443,12 @@ namespace MeshWorker matrix = &m; } + template + inline void + MGMatrixSimple::initialize(const MGConstraints& c) + { + mg_constraints = &c; + } template inline void @@ -1442,9 +1471,9 @@ namespace MeshWorker template - template + template inline void - MGMatrixSimple::initialize_info(DoFInfo& info, + MGMatrixSimple::initialize_info(DOFINFO& info, bool interior_face) const { info.initialize_matrices(1, interior_face); @@ -1458,17 +1487,75 @@ namespace MeshWorker MATRIX& G, const FullMatrix& M, const std::vector& i1, - const std::vector& i2) + const std::vector& i2, + const unsigned int level) { AssertDimension(M.m(), i1.size()); AssertDimension(M.n(), i2.size()); - for (unsigned int j=0; j= threshold) - G.add(i1[j], i2[k], M(j,k)); + if(mg_constraints == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(j,k)); + } + else + { + for (unsigned int j=0; j= threshold) + if(mg_constraints->at_refinement_edge(level, i1[j])==false && + mg_constraints->at_refinement_edge(level, i2[k])==false) + if((mg_constraints->is_boundary_index(level, i1[j])==false && + mg_constraints->is_boundary_index(level, i2[k])==false) + || + (mg_constraints->is_boundary_index(level, i1[j])==true && + mg_constraints->is_boundary_index(level, i2[k])==true + && + i1[j] == i2[k]) + ) + G.add(i1[j], i2[k], M(j,k)); + } } + template + inline void + MGMatrixSimple::assemble_down( + MATRIX& G, + const FullMatrix& M, + const std::vector& i1, + const std::vector& i2, + const unsigned int level) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + if(mg_constraints == 0) + { + for (unsigned int j=0; j= threshold) + G.add(i1[j], i2[k], M(j,k)); + } + else + { + for (unsigned int j=0; j= threshold) + if(mg_constraints->at_refinement_edge(level, i1[j])==true && + mg_constraints->at_refinement_edge(level, i2[k])==false) + if((mg_constraints->at_refinement_edge_boundary(level, i1[j])==false && + mg_constraints->at_refinement_edge_boundary(level, i2[k])==false) + || + (mg_constraints->at_refinement_edge_boundary(level, i1[j])==true && + mg_constraints->at_refinement_edge_boundary(level, i2[k])==true + && + i1[j] == i2[k]) + ) + G.add(i1[j], i2[k], M(j,k)); + } + } template inline void @@ -1476,33 +1563,59 @@ namespace MeshWorker MATRIX& G, const FullMatrix& M, const std::vector& i1, - const std::vector& i2) + const std::vector& i2, + const unsigned int level) { AssertDimension(M.n(), i1.size()); AssertDimension(M.m(), i2.size()); + if(mg_constraints == 0) + { for (unsigned int j=0; j= threshold) G.add(i1[j], i2[k], M(k,j)); + } + else + { + for (unsigned int j=0; j= threshold) + if(mg_constraints->at_refinement_edge(level, i1[j])==true && + mg_constraints->at_refinement_edge(level, i2[k])==false) + if((mg_constraints->at_refinement_edge_boundary(level, i1[j])==false && + mg_constraints->at_refinement_edge_boundary(level, i2[k])==false) + || + (mg_constraints->at_refinement_edge_boundary(level, i1[j])==true && + mg_constraints->at_refinement_edge_boundary(level, i2[k])==true + && + i1[j] == i2[k]) + ) + G.add(i1[j], i2[k], M(k,j)); + } } template - template + template inline void - MGMatrixSimple::assemble(const DoFInfo& info) + MGMatrixSimple::assemble(const DOFINFO& info) { const unsigned int level = info.cell->level(); - assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices); + assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); + if(mg_constraints != 0) + { + assemble_transpose((*interface_down)[level],info.matrix(0,false).matrix, info.indices, info.indices, level); + assemble_down((*interface_up)[level], info.matrix(0,false).matrix, info.indices, info.indices, level); + } } template - template + template inline void - MGMatrixSimple::assemble(const DoFInfo& info1, - const DoFInfo& info2) + MGMatrixSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) { const unsigned int level1 = info1.cell->level(); const unsigned int level2 = info2.cell->level(); @@ -1513,6 +1626,7 @@ namespace MeshWorker 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 { @@ -1520,9 +1634,9 @@ namespace MeshWorker // Do not add info2.M1, // which is done by // the coarser cell - 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); + 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); } } @@ -1561,10 +1675,10 @@ namespace MeshWorker template - template + template inline void MatrixLocalBlocksToGlobalBlocks::initialize_info( - DoFInfo& info, + DOFINFO& info, bool interior_face) const { info.initialize_matrices(*matrices, interior_face); @@ -1620,10 +1734,10 @@ namespace MeshWorker template - template + template inline void MatrixLocalBlocksToGlobalBlocks::assemble( - const DoFInfo& info) + const DOFINFO& info) { for (unsigned int i=0;isize();++i) { @@ -1638,11 +1752,11 @@ namespace MeshWorker template - template + template inline void MatrixLocalBlocksToGlobalBlocks::assemble( - const DoFInfo& info1, - const DoFInfo& info2) + const DOFINFO& info1, + const DOFINFO& info2) { for (unsigned int i=0;isize();++i) { @@ -1681,10 +1795,10 @@ namespace MeshWorker template - template + template inline void MGMatrixLocalBlocksToGlobalBlocks::initialize_info( - DoFInfo& info, + DOFINFO& info, bool interior_face) const { info.initialize_matrices(matrices, interior_face); @@ -1751,10 +1865,10 @@ namespace MeshWorker template - template + template inline void MGMatrixLocalBlocksToGlobalBlocks::assemble( - const DoFInfo& info) + const DOFINFO& info) { const unsigned int level = info.cell->level(); @@ -1772,11 +1886,11 @@ namespace MeshWorker template - template + template inline void MGMatrixLocalBlocksToGlobalBlocks::assemble( - const DoFInfo& info1, - const DoFInfo& info2) + const DOFINFO& info1, + const DOFINFO& info2) { const unsigned int level1 = info1.cell->level(); const unsigned int level2 = info2.cell->level(); @@ -1837,9 +1951,9 @@ namespace MeshWorker template - template + template inline void - SystemSimple::initialize_info(DoFInfo& info, + SystemSimple::initialize_info(DOFINFO& info, bool interior_face) const { MatrixSimple::initialize_info(info, interior_face); @@ -1848,9 +1962,9 @@ namespace MeshWorker template - template + template inline void - SystemSimple::assemble(const DoFInfo& info) + SystemSimple::assemble(const DOFINFO& info) { MatrixSimple::assemble(info); ResidualSimple::assemble(info); @@ -1858,10 +1972,10 @@ namespace MeshWorker template - template + template inline void - SystemSimple::assemble(const DoFInfo& info1, - const DoFInfo& info2) + SystemSimple::assemble(const DOFINFO& info1, + const DOFINFO& info2) { MatrixSimple::assemble(info1, info2); ResidualSimple::assemble(info1, info2); -- 2.39.5