symmetric DirichletBC and bug in FileWriter corrected

dune/amdis/DirichletBC.inc.hpp

@@ -31,11 +31,20 @@ namespace AMDiS
     using Dune::Functions::interpolate;
 
     test_exit_dbg(initialized, "Boundary condition not initialized!");
-    matrix.clearDirichletRows(dirichletNodes, apply);
+    auto columns = matrix.applyDirichletBC(dirichletNodes, apply);
+
 
     if (apply) {
       interpolate(matrix.getRowFeSpace(), wrapper(rhs.getVector()), values, dirichletNodes);
       interpolate(matrix.getColFeSpace(), wrapper(solution.getVector()), values, dirichletNodes);
+
+      // subtract columns of dirichlet nodes from rhs
+      for (size_t i = 0; i < columns.size(); ++i) {
+        for (auto const& aij : columns[i]) {
+          rhs[i] -= aij.second * solution[aij.first];
+        }
+      }
+
     }
   }
 

dune/amdis/FileWriter.hpp

@@ -128,7 +128,7 @@ namespace AMDiS
     }
 
   private:
-    MeshView const& meshView;
+    MeshView meshView;
     shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter;
     shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter;
 

dune/amdis/ProblemStat.hpp

@@ -114,6 +114,10 @@ namespace AMDiS
       : ProblemStatSeq(name)
     {
       this->mesh = std::shared_ptr<Mesh>(&mesh, optional_delete(false));
+      componentMeshes.resize(nComponents, this->mesh.get());
+
+      meshName = "";
+      Parameters::get(name + "->mesh", meshName);
     }
 
 
@@ -228,6 +232,7 @@ namespace AMDiS
     void setMesh(Mesh& mesh_)
     {
       mesh = std::shared_ptr<Mesh>(&mesh_, optional_delete(false));
+      std::fill(componentMeshes.begin(), componentMeshes.end(), this->mesh.get());
 
       createFeSpaces();
       createMatricesAndVectors();

dune/amdis/linear_algebra/istl/DOFMatrix.hpp

@@ -147,7 +147,7 @@ namespace AMDiS
       return matrix->nonzeroes();
     }
 
-    void clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
+    auto clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
     {
       // loop over the matrix rows
       for (std::size_t i = 0; i < matrix->N(); ++i) {
@@ -159,6 +159,9 @@ namespace AMDiS
                 *cIt = (apply && i == cIt.index() ? 1 : 0);
         }
       }
+
+      std::vector<std::map<size_t,value_type>> columns; // symmetric dbc not yet implemented
+      return columns;
     }
 
   private:

dune/amdis/linear_algebra/mtl/DOFMatrix.hpp

@@ -8,6 +8,8 @@
 #include <boost/numeric/mtl/utility/property_map.hpp>
 #include <boost/numeric/mtl/utility/range_wrapper.hpp>
 
+#include <dune/common/std/optional.hh>
+
 #include <dune/amdis/Output.hpp>
 #include <dune/amdis/common/ClonablePtr.hpp>
 
@@ -147,8 +149,10 @@ namespace AMDiS
 
     /// \brief Deletes all rows with \p dirichletNodes != 0 and if \p apply is true, adds
     /// a one on the diagonal of the matrix.
-    void clearDirichletRows(std::vector<char> const& dirichletNodes, bool apply)
+    auto applyDirichletBC(std::vector<char> const& dirichletNodes, bool apply)
     {
+      std::vector<std::map<size_t,value_type>> columns(num_rows(*matrix));
+
       // Define the property maps
       auto row   = mtl::mat::row_map(*matrix);
       auto col   = mtl::mat::col_map(*matrix);
@@ -157,14 +161,77 @@ namespace AMDiS
       // iterate over the matrix
       for (auto r : mtl::major_of(*matrix)) {   // rows or columns
         for (auto i : mtl::nz_of(r)) {          // non-zeros within
-          if (!dirichletNodes[row(i)])
-            break;
+          if (dirichletNodes[row(i)]) {
+            // set identity row
+            value(i, (apply && row(i) == col(i) ? value_type(1) : value_type(0)) );
+          } else if (apply && dirichletNodes[col(i)]) {
+            columns[row(i)][col(i)] = value(i);
+            value(i, value_type(0));
+          }
+        }
+      }
+
+      return columns;
+    }
+
+#if 0
+    void applyPeriodicBC(std::vector<char> const& periodicNodes,
+                         std::map<size_t, size_t> const& association, bool applyRow, bool applyCol)
+    {
+      bool apply = applyRow && applyCol;
+
+      // Define the property maps
+      auto row   = mtl::mat::row_map(*matrix);
+      auto col   = mtl::mat::col_map(*matrix);
+      auto value = mtl::mat::value_map(*matrix);
+
+      std::vector<std::map<size_t, std::list<value_type>>> row_values(num_cols(*matrix));
+      std::vector<std::map<size_t, std::list<value_type>>> col_values(num_rows(*matrix));
+      std::vector<char> dualNodes(periodicNodes.size(), 0);
+
+      // iterate over the matrix to collect the values and erase the source-row/col
+      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
+        for (auto i : mtl::nz_of(r)) {          // non-zeros within
+          if (applyRow && periodicNodes[row(i)]) {
+            row_values[col(i)][association[row(i)]].push_back(value(i));
+            value(i, value_type(0));
+            dualNodes[association[row(i)]] = 1;
+          } else if (applyCol && periodicNodes[col(i)]) {
+            col_values[row(i)][association[col(i)]].push_back(value(i));
+            value(i, value_type(0));
+          }
+        }
+      }
+
+      // TODO: use inserter for assignment of values!!!
+      // iterate over the matrix to assign the values
+      for (auto r : mtl::major_of(*matrix)) {   // rows or columns
+        for (auto i : mtl::nz_of(r)) {          // non-zeros within
+          if (applyRow && dualNodes[row(i)]) {
+            value(i, std::accumulate(row_values[col(i)][row(i)].begin(),
+                                     row_values[col(i)][row(i)].end(),
+                                     value(i)) );
+          }
+          if (applyCol && dualNodes[col(i)]) {
+            value(i, std::accumulate(col_values[row(i)][col(i)].begin(),
+                                     col_values[row(i)][col(i)].end(),
+                                     value(i)) );
+          }
+        }
+      }
 
-          // set identity row
-          value(i, (apply && row(i) == col(i) ? value_type(1) : value_type(0)) );
+      // assign values 1, -1 to diagonal and associated entries
+      if (apply) {
+        Inserter ins(*matrix);
+        for (auto const& a : association) {
+          ins[a.first][a.first]   << value_type( 1);
+          ins[a.second][a.second] << value_type( 1);
+          ins[a.first][a.second]  << value_type(-1);
+          ins[a.second][a.first]  << value_type(-1);
         }
       }
     }
+#endif
 
   protected:
     // Estimates the slot-size used in the inserter to reserve enough space per row.

src/gridview.cc

@@ -32,10 +32,6 @@ int main(int argc, char** argv)
   Dune::FieldVector<double, 2> upper = {{1.0, 1.0}};
   auto grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(lower, upper, n);
 
-  int ref = 0;
-  Parameters::get("elliptMesh->global refinements", ref);
-  grid->globalRefine(ref);
-
 
   ElliptProblem prob("ellipt", *grid);
   prob.initialize(INIT_ALL);