From: Joscha Gedicke <joscha.gedicke@iwr.uni-heidelberg.de>
Date: Sun, 6 Mar 2016 14:47:28 +0000 (+0100)
Subject: support for inhomogeneous constraints in MeshWorker::Assembler::SystemSimple
X-Git-Tag: v8.5.0-rc1~1212^2~7
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5a20033cc75c33ed73a9546fbae5b1d06e28abeb;p=dealii.git

support for inhomogeneous constraints in MeshWorker::Assembler::SystemSimple
---

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<class DOFINFO>
       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<class DOFINFO>
       void assemble(const DOFINFO &info1,
                     const DOFINFO &info2);
+    protected:
+      /**
+       * The vector of global matrices being assembled.
+       */
+      std::vector<SmartPointer<MatrixType,MatrixSimple<MatrixType> > > 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 <code>M</code> into the element at
@@ -213,22 +227,11 @@ namespace MeshWorker
                     const std::vector<types::global_dof_index> &i1,
                     const std::vector<types::global_dof_index> &i2);
 
-      /**
-       * The vector of global matrices being assembled.
-       */
-      std::vector<SmartPointer<MatrixType,MatrixSimple<MatrixType> > > matrix;
       /**
        * A pointer to the object containing constraints.
        */
       SmartPointer<const ConstraintMatrix,MatrixSimple<MatrixType> > 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<class DOFINFO>
       void assemble(const DOFINFO &info1,
                     const DOFINFO &info2);
+
+    private:
+      /**
+       * A pointer to the object containing constraints.
+       */
+      SmartPointer<const ConstraintMatrix,MatrixSimple<MatrixType> > constraints;
+      /**
+        * Assemble a single matrix <code>M</code> into the element at
+        * <code>index</code> in the vector #matrix.
+        */
+      void assemble(const FullMatrix<double> &M,
+                    const Vector<double> &vector,
+                    const unsigned int index,
+                    const std::vector<types::global_dof_index> &indices);
+
+      void assemble(const FullMatrix<double> &M,
+                    const Vector<double> &vector,
+                    const unsigned int index,
+                    const std::vector<types::global_dof_index> &i1,
+                    const std::vector<types::global_dof_index> &i2);
     };
 
 
@@ -1120,8 +1143,7 @@ namespace MeshWorker
     inline void
     SystemSimple<MatrixType,VectorType>::initialize(const ConstraintMatrix &c)
     {
-      MatrixSimple<MatrixType>::initialize(c);
-      ResidualSimple<VectorType>::initialize(c);
+      constraints = &c;
     }
 
 
@@ -1135,14 +1157,101 @@ namespace MeshWorker
       ResidualSimple<VectorType>::initialize_info(info, face);
     }
 
+    template <typename MatrixType,typename VectorType>
+    inline void
+    SystemSimple<MatrixType,VectorType>::assemble(const FullMatrix<double> &M,
+                                                  const Vector<double> &vector,
+                                                  const unsigned int index,
+                                                  const std::vector<types::global_dof_index> &indices)
+    {
+      AssertDimension(M.m(), indices.size());
+      AssertDimension(M.n(), indices.size());
+
+      AnyData residuals = ResidualSimple<VectorType>::residuals;
+      VectorType *v = residuals.entry<VectorType *>(index);
+
+      if (constraints == 0)
+        {
+          for (unsigned int i=0; i<indices.size(); ++i)
+            (*v)(indices[i]) += vector(i);
+
+          for (unsigned int j=0; j<indices.size(); ++j)
+            for (unsigned int k=0; k<indices.size(); ++k)
+              if (std::fabs(M(j,k)) >= MatrixSimple<MatrixType>::threshold)
+                MatrixSimple<MatrixType>::matrix[index]->add(indices[j], indices[k], M(j,k));
+        }
+      else
+        {
+          constraints->distribute_local_to_global(M,vector,indices,*MatrixSimple<MatrixType>::matrix[index],*v, true);
+        }
+    }
+
+    template <typename MatrixType,typename VectorType>
+    inline void
+    SystemSimple<MatrixType,VectorType>::assemble(const FullMatrix<double> &M,
+                                                  const Vector<double> &vector,
+                                                  const unsigned int index,
+                                                  const std::vector<types::global_dof_index> &i1,
+                                                  const std::vector<types::global_dof_index> &i2)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+
+      AnyData residuals = ResidualSimple<VectorType>::residuals;
+      VectorType *v = residuals.entry<VectorType *>(index);
+
+      if (constraints == 0)
+        {
+          for (unsigned int i=0; i<i1.size(); ++i)
+            (*v)(i1[i]) += vector(i);
+
+          for (unsigned int j=0; j<i1.size(); ++j)
+            for (unsigned int k=0; k<i2.size(); ++k)
+              if (std::fabs(M(j,k)) >= MatrixSimple<MatrixType>::threshold)
+                MatrixSimple<MatrixType>::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<MatrixType>::matrix[index]);
+        }
+    }
+
 
     template <typename MatrixType, typename VectorType>
     template <class DOFINFO>
     inline void
     SystemSimple<MatrixType,VectorType>::assemble(const DOFINFO &info)
     {
-      MatrixSimple<MatrixType>::assemble(info);
-      ResidualSimple<VectorType>::assemble(info);
+      AssertDimension(MatrixSimple<MatrixType>::matrix.size(),ResidualSimple<VectorType>::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<MatrixSimple<MatrixType>::matrix.size(); ++m)
+            assemble(info.matrix(m,false).matrix,info.vector(m).block(0), m, info.indices);
+        }
+      else
+        {
+          for (unsigned int m=0; m<MatrixSimple<MatrixType>::matrix.size(); ++m)
+            for (unsigned int k=0; k<n*n; ++k)
+              {
+                const unsigned int row = info.matrix(k+m*n*n,false).row;
+                const unsigned int column = info.matrix(k+m*n*n,false).column;
+
+                if (row == column)
+                  assemble(info.matrix(k+m*n*n,false).matrix,
+                           info.vector(m).block(column), m,
+                           info.indices_by_block[column]);
+                else
+                  assemble(info.matrix(k+m*n*n,false).matrix,
+                           info.vector(m).block(column), m,
+                           info.indices_by_block[row],
+                           info.indices_by_block[column]);
+              }
+        }
+
     }
 
 
@@ -1152,8 +1261,50 @@ namespace MeshWorker
     SystemSimple<MatrixType,VectorType>::assemble(const DOFINFO &info1,
                                                   const DOFINFO &info2)
     {
-      MatrixSimple<MatrixType>::assemble(info1, info2);
-      ResidualSimple<VectorType>::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<MatrixSimple<MatrixType>::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<MatrixSimple<MatrixType>::matrix.size(); ++m)
+            for (unsigned int k=0; k<n*n; ++k)
+              {
+                const unsigned int row = info1.matrix(k+m*n*n,false).row;
+                const unsigned int column = info1.matrix(k+m*n*n,false).column;
+
+                if (row == column)
+                  {
+                    assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(column), m,
+                             info1.indices_by_block[column]);
+                    assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(column), m,
+                             info2.indices_by_block[column]);
+                  }
+                else
+                  {
+                    assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(column), m,
+                             info1.indices_by_block[row], info1.indices_by_block[column]);
+                    assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(column), m,
+                             info2.indices_by_block[row], info2.indices_by_block[column]);
+                  }
+                assemble(info1.matrix(k+m*n*n,true).matrix, info1.vector(m).block(column), m,
+                         info1.indices_by_block[row], info2.indices_by_block[column]);
+                assemble(info2.matrix(k+m*n*n,true).matrix, info2.vector(m).block(column), m,
+                         info2.indices_by_block[row], info1.indices_by_block[column]);
+              }
+        }
     }
   }
 }