Improvement of implementation of improveGrid() method: Tetraeder mid points...

Improvement of implementation of improveGrid() method: Tetraeder mid points are moved, too if necessary; all grid levels are considered; grid improvement is iterated until no vertices are to be moved.

dune/fufem/improvegrid.hh

@@ -63,36 +63,21 @@ double aspectRatio(const EntityType& entity)
 }
 
 template <class GridType>
-void improveGrid(GridType& grid, double threshold)
+int improveGridLevel(GridType& grid, double threshold, int level)
 {
     const int dim = GridType::dimension;
-    int maxLevel = grid.maxLevel();
-
-    if (maxLevel==0)
-        return;
-
-    // Only works in 3d
-    if (dim!=3)
-        DUNE_THROW(Dune::GridError, "Grid improvement only for 3d grids!");
-
-    const auto& maxLevelView = grid.levelGridView(maxLevel);
-    const auto& indexSet = maxLevelView.indexSet();
+    const auto& levelView = grid.levelGridView(level);
+    const auto& indexSet = levelView.indexSet();
 
     Dune::BitSetVector<1> isNewOnThisLevel(indexSet.size(dim), false);
-    std::vector<Dune::FieldVector<double, dim> > vertexPosition(indexSet.size(dim));
+    Dune::BitSetVector<1> isTetMidPoint(indexSet.size(dim), false);
     std::vector<Dune::FieldVector<double, dim> > bisectionPosition(indexSet.size(dim));
 
-    // Copy all vertex positions into an indexable data structure
-    for (const auto& v : vertices(maxLevelView))
-        vertexPosition[indexSet.index(v)] = v.geometry().corner(0);
-
     // ///////////////////////////////////////////////////////////////
     //   Create the whole grid boundary as a boundary patch
     // ///////////////////////////////////////////////////////////////
 
-    BoundaryPatch<typename GridType::LevelGridView> boundary(maxLevelView, true);
-
-    std::cout << "Boundary contains " << boundary.numVertices() << " vertices\n";
+    BoundaryPatch<typename GridType::LevelGridView> boundary(levelView, true);
 
     // /////////////////////////////////////////////////////////////
     //   We need to know which vertices is new on the maxlevel and which
@@ -106,7 +91,7 @@ void improveGrid(GridType& grid, double threshold)
                                   {0.5,0,0.5},
                                   {0,0.5,0.5}};
 
-    for (const auto& e : elements(maxLevelView)) {
+    for (const auto& e : elements(levelView)) {
 
         const auto& father = e.father();
 
@@ -125,6 +110,15 @@ void improveGrid(GridType& grid, double threshold)
                     break;
                 }
             }
+            // tetraeder mid points must be moved, too
+            if ( std::abs(positionInFather[0]-0.25) < 0.1
+                 && std::abs(positionInFather[1]-0.25) < 0.1
+                 && std::abs(positionInFather[2]-0.25) < 0.1) {
+
+                isNewOnThisLevel[indexSet.subIndex(e,i,dim)][0] = true;
+                isTetMidPoint[indexSet.subIndex(e,i,dim)][0] = true;
+                bisectionPosition[indexSet.subIndex(e,i,dim)] = father.geometry().global(positionInFather);
+            }
         }
     }
 
@@ -132,7 +126,7 @@ void improveGrid(GridType& grid, double threshold)
     //   Mark the vertices of all 'bad' elements
     // ///////////////////////////////////////////////////////////////
     Dune::BitSetVector<1> badVertices(indexSet.size(dim),false);
-    for (const auto& e : elements(maxLevelView)){
+    for (const auto& e : elements(levelView)){
 
         Dune::FieldVector<double,dim> dummy(0.5);
 
@@ -141,17 +135,15 @@ void improveGrid(GridType& grid, double threshold)
 
             // This element is of bad quality.  Mark all its vertices
             for (size_t i=0; i<e.subEntities(dim); i++){
-                if (boundary.containsVertex(indexSet.subIndex(e,i,dim))) {
+                if (boundary.containsVertex(indexSet.subIndex(e,i,dim)) || isTetMidPoint[indexSet.subIndex(e,i,dim)][0]) {
                     badVertices[indexSet.subIndex(e,i,dim)][0] = true;
                 }
             }
         }
     }
 
-    std::cout << "The grid has " << badVertices.count() << " bad vertices!" << std::endl;
-
     // Loop over all new vertices
-    for (const auto& v : vertices(maxLevelView)) {
+    for (const auto& v : vertices(levelView)) {
 
         // Get the index
         int vIdx = indexSet.index(v);
@@ -161,6 +153,34 @@ void improveGrid(GridType& grid, double threshold)
 
         grid.setPosition(v, bisectionPosition[vIdx]);
     }
+
+    return badVertices.count();
+}
+
+template <class GridType>
+void improveGrid(GridType& grid, double threshold, int maxIter=10)
+{
+  const int dim = GridType::dimension;
+  int maxLevel = grid.maxLevel();
+
+  if (maxLevel==0)
+      return;
+
+  // Only works in 3d
+  if (dim!=3)
+      DUNE_THROW(Dune::GridError, "Grid improvement only for 3d grids!");
+
+  // By moving some vertices to improve certain elements other elements can worsen.
+  for(int i=0;i<maxIter;++i) {
+
+    int nBadVertices = 0;
+    // Could be, that new nodes were added on some corser level, too. So we must check all levels.
+    for(int level=maxLevel; level>0; --level) {
+      nBadVertices += improveGridLevel(grid, threshold, level);
+    }
+    std::cout << nBadVertices << " bad vertices were detected in grid improvement iteration " << i << "." << std::endl;
+    if(nBadVertices==0) break;
+  }
 }