Commit 22b5aba4 authored by Praetorius, Simon's avatar Praetorius, Simon

calcSignedDist for periodic boundary conditions

parent d6bd7df2
......@@ -21,6 +21,7 @@
#include "HL_SignedDistTraverse.h"
#include "VelocityExtFromVelocityField.h"
#include "VertexVector.h"
namespace reinit {
......@@ -146,6 +147,11 @@ void HL_SignedDistTraverse::HL_updateIteration()
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;
ElInfo *elInfo;
......@@ -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);
}
......
......@@ -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
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment