calcSignedDist for periodic boundary conditions

AMDiS/src/reinit/HL_SignedDistTraverse.cc

@@ -21,6 +21,7 @@
 
 #include "HL_SignedDistTraverse.h"
 #include "VelocityExtFromVelocityField.h"
+#include "VertexVector.h"
 
 namespace reinit {
 
@@ -145,6 +146,11 @@ void HL_SignedDistTraverse::HL_updateIteration()
     update_DOF = sD_DOF;
   else
     update_DOF = sDOld_DOF;
+  
+  std::vector<VertexVector*> associations;
+  for (int b : boundaryTypes) {
+    associations.push_back( &(feSpace->getMesh()->getPeriodicAssociations(BoundaryType{b}, feSpace->getAdmin())) );
+  }
 
   // ===== Iteration loop: proceed until tolerance is reached. =====
   TraverseStack stack;
@@ -159,8 +165,38 @@ void HL_SignedDistTraverse::HL_updateIteration()
     elInfo =
       stack.traverseFirst(feSpace->getMesh(), -1,
 			  Mesh::CALL_LEAF_EL | Mesh::FILL_BOUND | Mesh::FILL_COORDS);
+      
+    int dim = feSpace->getMesh()->getDim();
+    GeoIndex sideGeoIndex = INDEX_OF_DIM(dim - 1, dim);
+    std::vector<DegreeOfFreedom> localIndices(feSpace->getBasisFcts()->getNumber());
+    
     while (elInfo) {
       HL_elementUpdate(elInfo);
+      
+      // synch parallel DOFs
+      if (!boundaryTypes.empty())
+        feSpace->getBasisFcts()->getLocalIndices(elInfo->getElement(), feSpace->getAdmin(), localIndices);
+      
+      for (std::size_t i = 0; i < boundaryTypes.size(); ++i) {
+        BoundaryType boundaryType = boundaryTypes[i];
+        VertexVector& association = *associations[i];
+        
+        auto absmin = [](double v0, double v1) { return std::abs(v0) < std::abs(v1) ? v0 : v1; };
+        
+        for (int side = 0; side <= dim; side++) {
+          if (elInfo->getBoundary(sideGeoIndex, side) == boundaryType) {
+            for (int vertex = 0; vertex < dim; vertex++) {
+              DegreeOfFreedom right = localIndices[vertex];
+              DegreeOfFreedom left = association[right]; 
+          
+              
+              double value = absmin((*sD_DOF)[left], (*sD_DOF)[right]);
+              (*sD_DOF)[left] = (*sD_DOF)[right] = value;
+            }
+          }
+        }
+      }
+      
       elInfo = stack.traverseNext(elInfo);
     }
 

AMDiS/src/reinit/HL_SignedDistTraverse.h

@@ -86,6 +86,7 @@ public:
     Parameters::get(name + "->tolerance", tol);
     Parameters::get(name + "->maximal number of iteration steps", maxIt);
     Parameters::get(name + "->Gauss-Seidel iteration", GaussSeidelFlag);
+    Parameters::get(name + "->periodic boundaries", boundaryTypes);
 
     TEST_EXIT(tol > 0)("illegal tolerance !\n");
 
@@ -190,6 +191,8 @@ public:
   int setUpdate_Cntr;
   // ---> end: for test purposes
 
+  // boundary types for periodic mesh boundaries
+  std::vector<signed int> boundaryTypes;
 };
 
 } // end namespace reinit