From 5a20033cc75c33ed73a9546fbae5b1d06e28abeb Mon Sep 17 00:00:00 2001 From: Joscha Gedicke Date: Sun, 6 Mar 2016 15:47:28 +0100 Subject: [PATCH] support for inhomogeneous constraints in MeshWorker::Assembler::SystemSimple --- include/deal.II/meshworker/simple.h | 187 +++++++++++++++++++++++++--- 1 file changed, 169 insertions(+), 18 deletions(-) diff --git a/include/deal.II/meshworker/simple.h b/include/deal.II/meshworker/simple.h index bf8993b31c..9ea2e007f6 100644 --- a/include/deal.II/meshworker/simple.h +++ b/include/deal.II/meshworker/simple.h @@ -107,11 +107,12 @@ namespace MeshWorker template void assemble(const DOFINFO &info1, const DOFINFO &info2); - private: + protected: /** * The global residal vectors filled by assemble(). */ AnyData residuals; + private: /** * A pointer to the object containing constraints. */ @@ -203,6 +204,19 @@ namespace MeshWorker template void assemble(const DOFINFO &info1, const DOFINFO &info2); + protected: + /** + * The vector of global matrices being assembled. + */ + std::vector > > matrix; + + /** + * The smallest positive number that will be entered into the global + * matrix. All smaller absolute values will be treated as zero and will + * not be assembled. + */ + const double threshold; + private: /** * Assemble a single matrix M into the element at @@ -213,22 +227,11 @@ namespace MeshWorker const std::vector &i1, const std::vector &i2); - /** - * The vector of global matrices being assembled. - */ - std::vector > > matrix; /** * A pointer to the object containing constraints. */ SmartPointer > constraints; - /** - * The smallest positive number that will be entered into the global - * matrix. All smaller absolute values will be treated as zero and will - * not be assembled. - */ - const double threshold; - }; @@ -457,6 +460,26 @@ namespace MeshWorker template void assemble(const DOFINFO &info1, const DOFINFO &info2); + + private: + /** + * A pointer to the object containing constraints. + */ + SmartPointer > constraints; + /** + * Assemble a single matrix M into the element at + * index in the vector #matrix. + */ + void assemble(const FullMatrix &M, + const Vector &vector, + const unsigned int index, + const std::vector &indices); + + void assemble(const FullMatrix &M, + const Vector &vector, + const unsigned int index, + const std::vector &i1, + const std::vector &i2); }; @@ -1120,8 +1143,7 @@ namespace MeshWorker inline void SystemSimple::initialize(const ConstraintMatrix &c) { - MatrixSimple::initialize(c); - ResidualSimple::initialize(c); + constraints = &c; } @@ -1135,14 +1157,101 @@ namespace MeshWorker ResidualSimple::initialize_info(info, face); } + template + inline void + SystemSimple::assemble(const FullMatrix &M, + const Vector &vector, + const unsigned int index, + const std::vector &indices) + { + AssertDimension(M.m(), indices.size()); + AssertDimension(M.n(), indices.size()); + + AnyData residuals = ResidualSimple::residuals; + VectorType *v = residuals.entry(index); + + if (constraints == 0) + { + for (unsigned int i=0; i= MatrixSimple::threshold) + MatrixSimple::matrix[index]->add(indices[j], indices[k], M(j,k)); + } + else + { + constraints->distribute_local_to_global(M,vector,indices,*MatrixSimple::matrix[index],*v, true); + } + } + + template + inline void + SystemSimple::assemble(const FullMatrix &M, + const Vector &vector, + const unsigned int index, + const std::vector &i1, + const std::vector &i2) + { + AssertDimension(M.m(), i1.size()); + AssertDimension(M.n(), i2.size()); + + AnyData residuals = ResidualSimple::residuals; + VectorType *v = residuals.entry(index); + + if (constraints == 0) + { + for (unsigned int i=0; i= MatrixSimple::threshold) + MatrixSimple::matrix[index]->add(i1[j], i2[k], M(j,k)); + } + else + { + constraints->distribute_local_to_global(vector, i1, i2, *v, M); + constraints->distribute_local_to_global(M, i1, i2, *MatrixSimple::matrix[index]); + } + } + template template inline void SystemSimple::assemble(const DOFINFO &info) { - MatrixSimple::assemble(info); - ResidualSimple::assemble(info); + AssertDimension(MatrixSimple::matrix.size(),ResidualSimple::residuals.size()); + Assert(!info.level_cell, ExcMessage("Cell may not access level dofs")); + const unsigned int n = info.indices_by_block.size(); + + if (n == 0) + { + for (unsigned int m=0; m::matrix.size(); ++m) + assemble(info.matrix(m,false).matrix,info.vector(m).block(0), m, info.indices); + } + else + { + for (unsigned int m=0; m::matrix.size(); ++m) + for (unsigned int k=0; k::assemble(const DOFINFO &info1, const DOFINFO &info2) { - MatrixSimple::assemble(info1, info2); - ResidualSimple::assemble(info1, info2); + Assert(!info1.level_cell, ExcMessage("Cell may not access level dofs")); + Assert(!info2.level_cell, ExcMessage("Cell may not access level dofs")); + AssertDimension(info1.indices_by_block.size(),info2.indices_by_block.size()); + + const unsigned int n = info1.indices_by_block.size(); + + if (n == 0) + { + for (unsigned int m=0; m::matrix.size(); ++m) + { + assemble(info1.matrix(m,false).matrix, info1.vector(m).block(0), m, info1.indices); + assemble(info1.matrix(m,true).matrix, info1.vector(m).block(0), m, info1.indices, info2.indices); + assemble(info2.matrix(m,false).matrix, info2.vector(m).block(0), m, info2.indices); + assemble(info2.matrix(m,true).matrix, info2.vector(m).block(0), m, info2.indices, info1.indices); + } + } + else + { + for (unsigned int m=0; m::matrix.size(); ++m) + for (unsigned int k=0; k