From 2f1b1d6d9befa1fff21c79602e1573e820df098b Mon Sep 17 00:00:00 2001
From: Thomas Witkowski <thomas.witkowski@gmx.de>
Date: Thu, 17 Dec 2009 14:40:19 +0000
Subject: [PATCH] Added reinit library to AMDiS repository.
---
AMDiS/Reinit/Makefile | 44 +
AMDiS/Reinit/src/BoundaryElementDist.cc | 59 +
AMDiS/Reinit/src/BoundaryElementDist.h | 82 +
AMDiS/Reinit/src/BoundaryElementEdgeDist.cc | 55 +
AMDiS/Reinit/src/BoundaryElementEdgeDist.h | 37 +
.../Reinit/src/BoundaryElementLevelSetDist.cc | 21 +
.../Reinit/src/BoundaryElementLevelSetDist.h | 37 +
AMDiS/Reinit/src/BoundaryElementNormalDist.cc | 54 +
AMDiS/Reinit/src/BoundaryElementNormalDist.h | 47 +
AMDiS/Reinit/src/BoundaryElementTopDist.cc | 403 +++++
AMDiS/Reinit/src/BoundaryElementTopDist.h | 83 +
AMDiS/Reinit/src/ElementLevelSet.cc | 423 +++++
AMDiS/Reinit/src/ElementLevelSet.h | 497 ++++++
AMDiS/Reinit/src/ElementUpdate.h | 43 +
AMDiS/Reinit/src/ElementUpdate_2d.cc | 52 +
AMDiS/Reinit/src/ElementUpdate_2d.h | 33 +
AMDiS/Reinit/src/ElementUpdate_3d.cc | 609 ++++++++
AMDiS/Reinit/src/ElementUpdate_3d.h | 117 ++
AMDiS/Reinit/src/HL_SignedDist.cc | 375 +++++
AMDiS/Reinit/src/HL_SignedDist.h | 407 +++++
AMDiS/Reinit/src/HL_SignedDistBornemann.h | 916 +++++++++++
AMDiS/Reinit/src/HL_SignedDistLevels.cc | 1360 +++++++++++++++++
AMDiS/Reinit/src/HL_SignedDistLevels.h | 234 +++
AMDiS/Reinit/src/HL_SignedDistTraverse.cc | 275 ++++
AMDiS/Reinit/src/HL_SignedDistTraverse.h | 157 ++
AMDiS/Reinit/src/NormEps.cc | 29 +
AMDiS/Reinit/src/NormEps.h | 46 +
AMDiS/Reinit/src/VelocityExt.cc | 260 ++++
AMDiS/Reinit/src/VelocityExt.h | 258 ++++
.../src/VelocityExtFromVelocityField.cc | 42 +
.../Reinit/src/VelocityExtFromVelocityField.h | 181 +++
31 files changed, 7236 insertions(+)
create mode 100644 AMDiS/Reinit/Makefile
create mode 100644 AMDiS/Reinit/src/BoundaryElementDist.cc
create mode 100644 AMDiS/Reinit/src/BoundaryElementDist.h
create mode 100644 AMDiS/Reinit/src/BoundaryElementEdgeDist.cc
create mode 100644 AMDiS/Reinit/src/BoundaryElementEdgeDist.h
create mode 100644 AMDiS/Reinit/src/BoundaryElementLevelSetDist.cc
create mode 100644 AMDiS/Reinit/src/BoundaryElementLevelSetDist.h
create mode 100644 AMDiS/Reinit/src/BoundaryElementNormalDist.cc
create mode 100644 AMDiS/Reinit/src/BoundaryElementNormalDist.h
create mode 100644 AMDiS/Reinit/src/BoundaryElementTopDist.cc
create mode 100644 AMDiS/Reinit/src/BoundaryElementTopDist.h
create mode 100644 AMDiS/Reinit/src/ElementLevelSet.cc
create mode 100644 AMDiS/Reinit/src/ElementLevelSet.h
create mode 100644 AMDiS/Reinit/src/ElementUpdate.h
create mode 100644 AMDiS/Reinit/src/ElementUpdate_2d.cc
create mode 100644 AMDiS/Reinit/src/ElementUpdate_2d.h
create mode 100644 AMDiS/Reinit/src/ElementUpdate_3d.cc
create mode 100644 AMDiS/Reinit/src/ElementUpdate_3d.h
create mode 100644 AMDiS/Reinit/src/HL_SignedDist.cc
create mode 100644 AMDiS/Reinit/src/HL_SignedDist.h
create mode 100644 AMDiS/Reinit/src/HL_SignedDistBornemann.h
create mode 100644 AMDiS/Reinit/src/HL_SignedDistLevels.cc
create mode 100644 AMDiS/Reinit/src/HL_SignedDistLevels.h
create mode 100644 AMDiS/Reinit/src/HL_SignedDistTraverse.cc
create mode 100644 AMDiS/Reinit/src/HL_SignedDistTraverse.h
create mode 100644 AMDiS/Reinit/src/NormEps.cc
create mode 100644 AMDiS/Reinit/src/NormEps.h
create mode 100644 AMDiS/Reinit/src/VelocityExt.cc
create mode 100644 AMDiS/Reinit/src/VelocityExt.h
create mode 100644 AMDiS/Reinit/src/VelocityExtFromVelocityField.cc
create mode 100644 AMDiS/Reinit/src/VelocityExtFromVelocityField.h
diff --git a/AMDiS/Reinit/Makefile b/AMDiS/Reinit/Makefile
new file mode 100644
index 00000000..2d71747b
--- /dev/null
+++ b/AMDiS/Reinit/Makefile
@@ -0,0 +1,44 @@
+SHELL = /bin/sh
+COMPILE = g++
+DEFS = -DPACKAGE=\"AMDiS\" -DVERSION=\"0.1\" -DHAVE_DLFCN_H=1
+DEBUG = 0
+AMDIS_DIR = ../
+
+INCLUDES = -I./src -I$(AMDIS_DIR)/src -I$(AMDIS_DIR)/compositeFEM/src -I$(AMDIS_DIR)/lib/boost_1_34_1 -I$(AMDIS_DIR)/lib/mtl4
+
+ifeq ($(strip $(DEBUG)), 0)
+ CPPFLAGS = -O2
+else
+ CPPFLAGS = -g -O0
+endif
+CPPFLAGS += -ftemplate-depth-30
+
+VPATH = .:$(AMDIS_DIR)/Reinit/src
+
+LIBOFILES = \
+BoundaryElementDist.o \
+BoundaryElementNormalDist.o \
+BoundaryElementLevelSetDist.o \
+BoundaryElementTopDist.o \
+BoundaryElementEdgeDist.o \
+ElementLevelSet.o \
+ElementUpdate_2d.o \
+ElementUpdate_3d.o \
+HL_SignedDist.o \
+HL_SignedDistTraverse.o \
+VelocityExt.o \
+VelocityExtFromVelocityField.o \
+NormEps.o
+
+all: libreinit.a
+
+# statische Bibliothek
+libreinit.a: $(LIBOFILES)
+ rm -f lib/libreinit.a
+ ar cq lib/libreinit.a $(LIBOFILES)
+
+.cc.o: $*.cc
+ $(COMPILE) $(DEFS) $(INCLUDES) $(CPPFLAGS) -c -o $*.o $<
+
+clean:
+ -rm -rf *.o
diff --git a/AMDiS/Reinit/src/BoundaryElementDist.cc b/AMDiS/Reinit/src/BoundaryElementDist.cc
new file mode 100644
index 00000000..b3a4adde
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementDist.cc
@@ -0,0 +1,59 @@
+#include "BoundaryElementDist.h"
+
+void
+BoundaryElementDist::calcNormal(
+ const FixVec<WorldVector<double>, DIMEN> &vecs,
+ WorldVector<double> &normalVec)
+{
+ FUNCNAME("BoundaryElementDist::calcNormal");
+
+ switch(dim) {
+ case 2: calcNormal_2d(vecs, normalVec);
+ break;
+ case 3: calcNormal_3d(vecs, normalVec);
+ break;
+ default: ERROR_EXIT("illegal dimension !\n");
+ }
+}
+
+void
+BoundaryElementDist::calcNormal_2d(
+ const FixVec<WorldVector<double>, DIMEN> &vecs,
+ WorldVector<double> &normalVec)
+{
+ FUNCNAME("BoundaryElementDist::calcNormal_2d");
+
+ TEST_EXIT(vecs.size() == dim)("illegal number of world vectors in vecs !\n");
+
+ double norm2 = 0.0;
+ double val;
+ for (int i=0; i<dim; ++i){
+ val = vecs[0][i] - vecs[1][i];
+ norm2 += val*val;
+ normalVec[dim-1-i] = val;
+ }
+ normalVec[0] *= -1;
+ double norm = sqrt(norm2);
+ for(int i=0; i<dim; ++i) {
+ normalVec[i] = 1/norm * normalVec[i];
+ }
+}
+
+void
+BoundaryElementDist::calcNormal_3d(
+ const FixVec<WorldVector<double>, DIMEN> &vecs,
+ WorldVector<double> &normalVec)
+{
+ FUNCNAME("BoundaryElementDist::calcNormal_3d");
+
+ TEST_EXIT(vecs.size() == dim)("illegal number of world vectors in vecs !\n");
+
+ WorldVector<double> A = vecs[1]-vecs[0];
+ WorldVector<double> B = vecs[2]-vecs[0];
+ vectorProduct(A, B, normalVec);
+
+ double norm = sqrt(normalVec * normalVec);
+ for(int i=0; i<dim; ++i) {
+ normalVec[i] = 1/norm * normalVec[i];
+ }
+}
diff --git a/AMDiS/Reinit/src/BoundaryElementDist.h b/AMDiS/Reinit/src/BoundaryElementDist.h
new file mode 100644
index 00000000..5088edfd
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementDist.h
@@ -0,0 +1,82 @@
+#ifndef BOUNDARYELEMENTDIST_H
+#define BOUNDARYELEMENTDIST_H
+
+#include "ElInfo.h"
+#include "FixVec.h"
+
+#include "ElementLevelSet.h"
+
+using namespace AMDiS;
+
+class BoundaryElementDist {
+ public:
+ MEMORY_MANAGED(BoundaryElementDist);
+
+ /**
+ * Constructor.
+ */
+ BoundaryElementDist(ElementLevelSet *elLS_, int dim_)
+ : dim(dim_),
+ elLS(elLS_)
+ {
+ FUNCNAME("BoundaryElementDist::BoundaryElementDist()");
+
+ TEST_EXIT(dim == 2 || dim == 3)
+ ("function only works for dimension 2 !\n");
+ };
+
+ /**
+ * Virtual destructor.
+ */
+ virtual ~BoundaryElementDist() {};
+
+ /**
+ * Calculates distance from the interface for all vertices of a boundary
+ * element.
+ *
+ * Pure virtual: Realizations in
+ * BoundaryElementLevelSetDist - boundary value initialization
+ * by level set function
+ * BoundaryElementTopDist - topological distance
+ * (distance within element)
+ * BoundaryElementEdgeDist - distance along edges
+ * BoundaryElementNormalDist - normal distance to
+ * intersection plane
+ */
+ virtual int calcDistOnBoundaryElement(ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec) = 0;
+
+ protected:
+ /**
+ * Calculate the (unit) normal for the plane defined by the world vectors in
+ * vecs.
+ */
+ void calcNormal(const FixVec<WorldVector<double>, DIMEN> &vecs,
+ WorldVector<double> &normalVec);
+
+ /**
+ * Normal calculation: 2-dimensional case.
+ */
+ void calcNormal_2d(const FixVec<WorldVector<double>, DIMEN> &vecs,
+ WorldVector<double> &normalVec);
+
+ /**
+ * Normal calculation: 3-dimensional case.
+ */
+ void calcNormal_3d(const FixVec<WorldVector<double>, DIMEN> &vecs,
+ WorldVector<double> &normalVec);
+
+ protected:
+ /**
+ * Dimension.
+ */
+ int dim;
+
+ /**
+ * Holds level set function and functionalities for intersection point
+ * calculation.
+ */
+ ElementLevelSet *elLS;
+};
+
+#endif // BOUNDARYELEMENTDIST_H
diff --git a/AMDiS/Reinit/src/BoundaryElementEdgeDist.cc b/AMDiS/Reinit/src/BoundaryElementEdgeDist.cc
new file mode 100644
index 00000000..73d76a5d
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementEdgeDist.cc
@@ -0,0 +1,55 @@
+#include "BoundaryElementEdgeDist.h"
+
+int
+BoundaryElementEdgeDist::calcDistOnBoundaryElement(
+ ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec)
+{
+ WorldVector<double> vert1, vert2;
+ int vert1Ind, vert2Ind;
+ double length;
+ double val;
+
+
+ // Get intersection information.
+ int elStatus = elLS->createElementLevelSet(elInfo);
+ if (elStatus != ElementLevelSet::LEVEL_SET_BOUNDARY)
+ return elStatus;
+
+ VectorOfFixVecs<DimVec<double> > *elIntersecPoints =
+ elLS->getElIntersecPoints();
+ int numIntersecPoints = elLS->getNumElIntersecPoints();
+
+ // For all vertices: calculate distance to intersection points.
+ for (int i = 0; i < numIntersecPoints; ++i) {
+ vert1Ind = -1;
+ vert2Ind = -1;
+ length = 0;
+
+ // Get vertices on edge with intersection point.
+ for (int j=0; j<=dim; ++j) {
+ if ((*elIntersecPoints)[i][j] > 1.e-15) {
+ if (vert1Ind == -1) vert1Ind = j;
+ else {
+ vert2Ind = j;
+ break;
+ }
+ }
+ }
+
+ // Get length of edge with intersection point.
+ vert1 = elInfo->getCoord(vert1Ind);
+ vert2 = elInfo->getCoord(vert2Ind);
+ for (int j=0; j<dim; ++j) {
+ length += (vert1[j] - vert2[j])*(vert1[j]-vert2[j]);
+ }
+ length = sqrt(length);
+
+ // Get distance of vert1 and vert2 to intersection point.
+ val = (*elIntersecPoints)[i][vert2Ind] * length;
+ dVec[vert1Ind] = val;
+ dVec[vert2Ind] = length - val;
+ }
+
+ return elStatus;
+}
diff --git a/AMDiS/Reinit/src/BoundaryElementEdgeDist.h b/AMDiS/Reinit/src/BoundaryElementEdgeDist.h
new file mode 100644
index 00000000..579fa634
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementEdgeDist.h
@@ -0,0 +1,37 @@
+#ifndef BOUNDARYELEMENTEDGEDIST_H
+#define BOUNDARYELEMENTEDGEDIST_H
+
+#include "ElInfo.h"
+#include "FixVec.h"
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+
+using namespace AMDiS;
+
+class BoundaryElementEdgeDist : public BoundaryElementDist
+{
+ public:
+ MEMORY_MANAGED(BoundaryElementEdgeDist);
+
+ /**
+ * Constructor.
+ */
+ BoundaryElementEdgeDist(ElementLevelSet *elLS_,
+ int dim_)
+ : BoundaryElementDist(elLS_, dim_)
+ {};
+
+ /**
+ * Calculates distance from the interface for all vertices of a boundary
+ * element.
+ * Distance is here the distance along edges.
+ *
+ * Return value: Status of element elInfo.
+ */
+ int calcDistOnBoundaryElement(ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec);
+};
+
+#endif // BOUNDARYELEMENTEDGEDIST_H
diff --git a/AMDiS/Reinit/src/BoundaryElementLevelSetDist.cc b/AMDiS/Reinit/src/BoundaryElementLevelSetDist.cc
new file mode 100644
index 00000000..69ad306e
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementLevelSetDist.cc
@@ -0,0 +1,21 @@
+#include "BoundaryElementLevelSetDist.h"
+
+int
+BoundaryElementLevelSetDist::calcDistOnBoundaryElement(
+ ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec)
+{
+ // Get intersection information.
+ int elStatus = elLS->createElementLevelSet(elInfo);
+ if (elStatus != ElementLevelSet::LEVEL_SET_BOUNDARY)
+ return elStatus;
+
+ const double *elVertLevelSetVec = elLS->getElVertLevelSetVec();
+
+ // Set distance to values of level set function in element vertices.
+ for (int i=0; i<=dim; ++i) {
+ dVec[i] = fabs(elVertLevelSetVec[i]);
+ }
+
+ return elStatus;
+}
diff --git a/AMDiS/Reinit/src/BoundaryElementLevelSetDist.h b/AMDiS/Reinit/src/BoundaryElementLevelSetDist.h
new file mode 100644
index 00000000..f639e207
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementLevelSetDist.h
@@ -0,0 +1,37 @@
+#ifndef BOUNDARYELEMENTLEVELSETDIST_H
+#define BOUNDARYELEMENTLEVELSETDIST_H
+
+#include "ElInfo.h"
+#include "FixVec.h"
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+
+using namespace AMDiS;
+
+class BoundaryElementLevelSetDist : public BoundaryElementDist
+{
+ public:
+ MEMORY_MANAGED(BoundaryElementLevelSetDist);
+
+ /**
+ * Constructor.
+ */
+ BoundaryElementLevelSetDist(ElementLevelSet *elLS_, int dim_)
+ : BoundaryElementDist(elLS_, dim_)
+ {};
+
+ /**
+ * Calculates distance from the interface for all vertices of a boundary
+ * element.
+ * Distance is here the value of the level set function at boundary
+ * vertices..
+ *
+ * Return value: Status of element elInfo.
+ */
+ int calcDistOnBoundaryElement(ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec);
+};
+
+#endif // BOUNDARYELEMENTLEVELSETDIST_H
diff --git a/AMDiS/Reinit/src/BoundaryElementNormalDist.cc b/AMDiS/Reinit/src/BoundaryElementNormalDist.cc
new file mode 100644
index 00000000..c14f5b93
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementNormalDist.cc
@@ -0,0 +1,54 @@
+#include "BoundaryElementNormalDist.h"
+
+int
+BoundaryElementNormalDist::calcDistOnBoundaryElement(
+ ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec)
+{
+
+ // Get intersection information.
+ int elStatus = elLS->createElementLevelSet(elInfo);
+ if (elStatus != ElementLevelSet::LEVEL_SET_BOUNDARY)
+ return elStatus;
+
+ VectorOfFixVecs<DimVec<double> > *elIntersecPoints =
+ elLS->getElIntersecPoints();
+
+ // Calculate (unit) normal to intersection plane.
+ FixVec<WorldVector<double>, DIMEN> planeVecs(dim, NO_INIT);
+ WorldVector<double> normalVec;
+ for (int i=0; i<dim; ++i) {
+ elInfo->coordToWorld((*elIntersecPoints)[i], (planeVecs[i]));
+ }
+ calcNormal(planeVecs, normalVec);
+
+ // Calculate normal distance for all vertices.
+ for (int i=0; i<=dim; ++i) {
+
+ dVec[i] = calculateDistLevelSetWithNormal(elInfo,
+ i,
+ normalVec,
+ planeVecs[0]);
+ }
+
+ return elStatus;
+}
+
+double
+BoundaryElementNormalDist::calculateDistLevelSetWithNormal(
+ ElInfo *elInfo,
+ int vert,
+ WorldVector<double> &normalVec,
+ WorldVector<double> &planeVec)
+{
+ // Get world coordinates of vertex.
+ const WorldVector<double> &vertex = elInfo->getCoord(vert);
+
+ // Calculate distance.
+ double dist = 0.0;
+ for (int i=0; i<dim; ++i) {
+ dist += (vertex[i]-planeVec[i]) * normalVec[i];
+ }
+
+ return fabs(dist);
+}
diff --git a/AMDiS/Reinit/src/BoundaryElementNormalDist.h b/AMDiS/Reinit/src/BoundaryElementNormalDist.h
new file mode 100644
index 00000000..a74c2725
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementNormalDist.h
@@ -0,0 +1,47 @@
+#ifndef BOUNDARYELEMENTNORMALDIST_H
+#define BOUNDARYELEMENTNORMALDIST_H
+
+#include "ElInfo.h"
+#include "FixVec.h"
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+
+using namespace AMDiS;
+
+class BoundaryElementNormalDist : public BoundaryElementDist
+{
+ public:
+ MEMORY_MANAGED(BoundaryElementNormalDist);
+
+ /**
+ * Constructor.
+ */
+ BoundaryElementNormalDist(ElementLevelSet *elLS_, int dim_)
+ : BoundaryElementDist(elLS_, dim_)
+ {};
+
+ /**
+ * Calculates distance from the interface for all vertices of a boundary
+ * element.
+ * Distance is here the normal distance.
+ *
+ * Return value: Status of element elInfo.
+ */
+ int calcDistOnBoundaryElement(ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec);
+
+ protected:
+ /**
+ * Calculates the distance of vertex vert in element elInfo
+ * from the plane given by the normal vector of the plane normalVec
+ * and one point lying in the plane planeVec.
+ */
+ double calculateDistLevelSetWithNormal(ElInfo *elInfo,
+ int vert,
+ WorldVector<double> &normalVec,
+ WorldVector<double> &planeVec);
+};
+
+#endif // BOUNDARYELEMENTNORMALDIST_H
diff --git a/AMDiS/Reinit/src/BoundaryElementTopDist.cc b/AMDiS/Reinit/src/BoundaryElementTopDist.cc
new file mode 100644
index 00000000..367337fc
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementTopDist.cc
@@ -0,0 +1,403 @@
+#include "BoundaryElementTopDist.h"
+
+int
+BoundaryElementTopDist::calcDistOnBoundaryElement(
+ ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec)
+{
+ //Cartesian coordinates of the 3 or 4 vertices of a simplex
+ FixVec<WorldVector<double>,VERTEX> Vert(dim, NO_INIT);
+ //normal unit vector in Cartesian coordinates
+ WorldVector<double> normalVec;
+ //Cartesian coordinates of dim intersection points
+ FixVec<WorldVector<double>,DIMEN> planeVecs_DIMEN(dim, NO_INIT);
+ //Cartesian coordinates of 3 or 4 intersection points (3d)
+ FixVec<WorldVector<double>,VERTEX> planeVecs(dim, NO_INIT);
+ //barycentric coordinates of the intersection point
+ //(straight line - normal)
+ DimVec<double> SP_Vec(dim, NO_INIT);
+ //Cartesian coordinates of SP_Vec
+ WorldVector<double> SP;
+
+ double lambda;
+ double h_dist = 0.0;
+ double h_dist2 = 0.0;
+ double dist = 0.0;
+
+ // for VelocityExt
+ int edgeUpdateInd = -1;
+ int sP1NextEdge = -1;
+ int sP2NextEdge = -1;
+ double lambdaNextPt = 0.0;
+
+ // Get intersection information.
+ int elStatus = elLS->createElementLevelSet(elInfo);
+ if (elStatus != ElementLevelSet::LEVEL_SET_BOUNDARY)
+ return elStatus;
+
+ //bar. coordinates of the intersection points
+ VectorOfFixVecs<DimVec<double> > *elIntersecPoints =
+ elLS->getElIntersecPoints();
+
+ //Cartesian coordinates of the vertices
+ for( int i=0; i<=dim; i++)
+ {
+ Vert[i] = elInfo->getCoord(i);
+ }
+
+ //Cartesian coordinates of the intersection point
+ //attention: in 3d it is possible that there are 4 intersection
+ //points, but her we save only 3 of them
+ for(int i=0; i<dim/*ElementLevelSet::getNumElIntersecPoints()*/; i++)
+ {
+ elInfo->coordToWorld((*elIntersecPoints)[i], (planeVecs_DIMEN[i]));
+ }
+
+ //calculating of the normal unit vector
+ calcNormal(planeVecs_DIMEN, normalVec);
+
+ //2D or 3D
+ if( dim == 2)
+ {
+ //loop over all vertices
+ for(int i=0; i<=dim; i++)
+ {
+ //Cartesian coordinates of the intersection point
+ //(straight line- normal)
+ h_dist = 0.0;
+ for (int j=0; j<dim; ++j)
+ {
+ h_dist += (Vert[i][j]-planeVecs_DIMEN[0][j]) * normalVec[j];
+ }
+ SP[0]=Vert[i][0]-h_dist*normalVec[0];
+ SP[1]=Vert[i][1]-h_dist*normalVec[1];
+
+ //barycentric coordinates of the intersection point
+ elInfo->worldToCoord(SP, &SP_Vec);
+
+ //calculating of the distance
+ //intersection point inside of the 2d-simplex?
+ if(SP_Vec[0]>=0 && SP_Vec[1]>=0 && SP_Vec[2]>=0)
+ {
+ dist = fabs(h_dist);
+ dVec[i] = dist;
+ //save barycentric coordinates of intersection point,
+ //for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_2D_boundary(SP_Vec[0], SP_Vec[1], SP_Vec[2], i);
+ }
+ }
+ //dist = min of the distance to the two intersection points
+ else
+ {
+ for (int j=0; j<elLS->getNumElIntersecPoints(); j++)
+ {
+ h_dist = 0.0;
+ for(int k=0; k<dim; k++)
+ {
+ h_dist += (Vert[i][k]-planeVecs_DIMEN[j][k])*
+ (Vert[i][k]-planeVecs_DIMEN[j][k]);
+ }
+ if (j == 0)
+ {
+ h_dist2 = h_dist;
+
+ // Store index of intersection point which is next
+ // point on interface.
+ edgeUpdateInd = j;
+ }
+ else if (h_dist < h_dist2)
+ {
+ h_dist2 = h_dist;
+
+ // Store index of intersection point which is next
+ // point on interface.
+ edgeUpdateInd = j;
+ }
+ }
+ dist = sqrt(h_dist2);
+ dVec[i] = dist;
+ //save barycentric coordinates of intersection point,
+ //for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_2D_boundary(
+ (*elIntersecPoints)[edgeUpdateInd][0],
+ (*elIntersecPoints)[edgeUpdateInd][1],
+ (*elIntersecPoints)[edgeUpdateInd][2],
+ i);
+ }
+ }//end of the else-part according to "if(!normalDistCalc)"
+ }//end of the loop over all vertices
+ }//end of the 2d-case
+
+ else //dim ==3
+ {
+ //Cartesian coordinates of the intersection point
+ //attention: in 3d it is possible that there are 4 intersection
+ //points, here we save them all
+ for(int i=0; i<elLS->getNumElIntersecPoints(); i++)
+ {
+ elInfo->coordToWorld((*elIntersecPoints)[i], (planeVecs[i]));
+ }
+
+ //loop over all 4 vertices
+ for ( int i=0; i<=dim; i++)
+ {
+ h_dist = 0.0;
+
+ //intersection point of the normal with the intersection plane
+ for ( int j = 0; j<dim; j++)
+ {
+ h_dist+=(Vert[i][j]-planeVecs[0][j])*normalVec[j];
+ }
+ for ( int j=0; j<dim; j++)
+ {
+ SP[j]=Vert[i][j]-h_dist*normalVec[j];
+ }
+
+ //barycentric coordinates of the intersection point
+ elInfo->worldToCoord(SP, &SP_Vec);
+
+ //intersection point inside of the 3d-simplex
+ //distance calculated with the unit normal vector of the
+ //intersection plane
+ if (SP_Vec[0]>=0 && SP_Vec[1]>=0 && SP_Vec[2]>=0 && SP_Vec[3]>=0)
+ {
+ dist = fabs(h_dist);
+ dVec[i] = dist;
+ //save barycentric coordinates of intersection point,
+ //for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_3D_boundary(SP_Vec[0], SP_Vec[1], SP_Vec[2], SP_Vec[3], i);
+ }
+ }
+ //dist = min of the distances to the cutting straight lines at
+ //the sides of the simplex
+ else if(elLS->getNumElIntersecPoints() == 3)
+ {
+ dist = 1.e15;
+ for (int j=0; j<dim; j++)
+ {
+ h_dist = calcDistOnIntersecEdge(planeVecs[j],
+ planeVecs[(j+1)%3],
+ Vert[i],
+ lambda);
+ // dist = (dist < h_dist) ? dist : h_dist;
+ if(h_dist<dist)
+ {
+ dist = h_dist;
+ if (velExt != NULL)
+ {
+ sP1NextEdge = j;
+ sP2NextEdge = (j+1)%3;
+ lambdaNextPt = lambda;
+ }
+ }
+ // if (velExt != NULL) {
+// sP1NextEdge = j;
+// sP2NextEdge = (j+1)%3;
+// lambdaNextPt = lambda;
+
+ }
+ dVec[i] = dist;
+ //save barycentric coordinates of intersection point,
+ //for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->calcBarycentricCoords_3D_boundary(
+ (*elIntersecPoints)[sP1NextEdge],
+ (*elIntersecPoints)[sP2NextEdge],
+ lambdaNextPt,
+ i);
+ }
+ }
+ else if(elLS->getNumElIntersecPoints() == 4)
+ {
+ dist = 1.e15;
+ for (int j=0; j<=1; j++)
+ {
+ h_dist = calcDistOnIntersecEdge(planeVecs[j],
+ planeVecs[(j+2)%3],
+ Vert[i],
+ lambda);
+ // dist = (dist < h_dist) ? dist : h_dist;
+
+// if (velExt != NULL) {
+// sP1NextEdge = j;
+// sP2NextEdge = (j+2)%3;
+// lambdaNextPt = lambda;
+// }
+ if(h_dist<dist)
+ {
+ dist = h_dist;
+ if (velExt != NULL)
+ {
+ sP1NextEdge = j;
+ sP2NextEdge = (j+2)%3;
+ lambdaNextPt = lambda;
+ }
+ }
+ }
+
+ h_dist = calcDistOnIntersecEdge(planeVecs[1],
+ planeVecs[2],
+ Vert[i],
+ lambda);
+ // dist = (dist < h_dist) ? dist : h_dist;
+
+// if (velExt != NULL) {
+// sP1NextEdge = 1;
+// sP2NextEdge = 2;
+// lambdaNextPt = lambda;
+// }
+ if(h_dist<dist)
+ {
+ dist = h_dist;
+ if (velExt != NULL)
+ {
+ sP1NextEdge = 1;
+ sP2NextEdge = 2;
+ lambdaNextPt = lambda;
+ }
+ }
+
+ h_dist = calcDistOnIntersecEdge(planeVecs[0],
+ planeVecs[3],
+ Vert[i],
+ lambda);
+ // dist = (dist < h_dist) ? dist : h_dist;
+
+// if (velExt != NULL) {
+// sP1NextEdge = 0;
+// sP2NextEdge = 3;
+// lambdaNextPt = lambda;
+// }
+ if(h_dist<dist)
+ {
+ dist = h_dist;
+ if (velExt != NULL)
+ {
+ sP1NextEdge = 0;
+ sP2NextEdge = 3;
+ lambdaNextPt = lambda;
+ }
+ }
+
+ dVec[i] = dist;
+ //save barycentric coordinates of intersection point,
+ //for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->calcBarycentricCoords_3D_boundary(
+ (*elIntersecPoints)[sP1NextEdge],
+ (*elIntersecPoints)[sP2NextEdge],
+ lambdaNextPt,
+ i);
+ }
+ }
+ }//end of the loop over all 4 vertices
+ }//end of the else-part according to "if( dim == 2)"
+
+ return elStatus;
+}
+
+double
+BoundaryElementTopDist::calcDistOnIntersecEdge(
+ const WorldVector<double> &sP1,
+ const WorldVector<double> &sP2,
+ const WorldVector<double> &v,
+ double &lambda)
+{
+ WorldVector<double> SP_projected;
+ double dist = 0.0;
+
+ // Calculate projection of v on straight line through sP1 and sP2.
+ projected_on_a_straight_line(sP1, sP2, v, SP_projected, lambda);
+
+ // Get distance and next point on interface (implicitly: lambda).
+ cases_for_lambda(v, SP_projected, sP1, sP2, lambda, dist);
+
+ return dist;
+}
+
+double
+BoundaryElementTopDist::calc_dist_between_two_points (
+ const WorldVector<double> &point1,
+ const WorldVector<double> &point2)
+{
+ double h_dist = 0.0;
+
+ for(int i=0; i<dim; i++)
+ {
+ h_dist += (point1[i] - point2[i])*(point1[i] - point2[i]);
+ }
+ h_dist = sqrt(h_dist);
+
+ return h_dist;
+}
+
+void
+BoundaryElementTopDist::cases_for_lambda(
+ const WorldVector<double> &SP,
+ const WorldVector<double> &SP_projected,
+ const WorldVector<double> &point1,
+ const WorldVector<double> &point2,
+ double &lambda,
+ double &dist)
+{
+ if(0 <= lambda && lambda <= 1)
+ {
+ dist = calc_dist_between_two_points(SP, SP_projected);
+ }
+ if(lambda < 0.0)
+ {
+ lambda = 0;
+ dist = calc_dist_between_two_points(SP, point1);
+ }
+ if(1.0 < lambda)
+ {
+ lambda = 1;
+ dist = calc_dist_between_two_points(SP, point2);
+ }
+}
+
+void
+BoundaryElementTopDist::projected_on_a_straight_line(
+ const WorldVector<double> &sP1,
+ const WorldVector<double> &sP2,
+ const WorldVector<double> &v,
+ WorldVector<double> &v_p,
+ double &lambda)
+ // lambda < 0: v_p---sP1---sP2
+ // 0<lambda<1: sP1---v_p---sP2
+ // 1 < lambda: sP1---sP2---v_p
+{
+ lambda = 0.0;
+ WorldVector<double> sP1_sP2;
+ WorldVector<double> sP1_v;
+ double norm_sP1_sP2 = 0.0;
+
+ for (int i=0; i<dim; ++i) {
+ sP1_sP2[i]=sP2[i]-sP1[i];
+ }
+ for (int i=0; i<dim; ++i) {
+ sP1_v[i]=v[i]-sP1[i];
+ }
+ for (int i=0; i<dim; ++i) {
+ norm_sP1_sP2 += sP1_sP2[i]*sP1_sP2[i];
+ }
+ // norm_sP1_sP2 = sqrt(norm_sP1_sP2);
+
+ for (int i=0; i<dim; ++i) {
+ lambda += sP1_v[i]*sP1_sP2[i];
+ }
+ lambda = lambda / norm_sP1_sP2;
+
+ for (int i=0; i<dim; i++)
+ {
+ v_p[i] = sP1[i] + lambda * sP1_sP2[i];
+ }
+}
diff --git a/AMDiS/Reinit/src/BoundaryElementTopDist.h b/AMDiS/Reinit/src/BoundaryElementTopDist.h
new file mode 100644
index 00000000..0adfb23b
--- /dev/null
+++ b/AMDiS/Reinit/src/BoundaryElementTopDist.h
@@ -0,0 +1,83 @@
+#ifndef BOUNDARYELEMENTTOPDIST_H
+#define BOUNDARYELEMENTTOPDIST_H
+
+#include "ElInfo.h"
+#include "FixVec.h"
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+#include "VelocityExt.h"
+
+using namespace AMDiS;
+
+class BoundaryElementTopDist : public BoundaryElementDist
+{
+ public:
+ MEMORY_MANAGED(BoundaryElementTopDist);
+
+ /**
+ * Constructor.
+ */
+ BoundaryElementTopDist(ElementLevelSet *elLS_,
+ int dim_,
+ VelocityExt *velExt_ = NULL)
+ : BoundaryElementDist(elLS_, dim_),
+ velExt(velExt_)
+ {};
+
+ /**
+ * Destructor.
+ */
+ ~BoundaryElementTopDist()
+ {};
+
+ /**
+ * Calculates distance from the interface for all vertices of a boundary
+ * element.
+ * Distance is here the topological distance.
+ *
+ * Return value: Status of element elInfo.
+ */
+ int calcDistOnBoundaryElement(ElInfo *elInfo,
+ FixVec<double, VERTEX> &dVec);
+
+ protected:
+ /**
+ * Calculates distance of vertex v to intersection edge through
+ * sP1 and sP2 within element.
+ * lambda implicitly gives the next point on the intersection edge.
+ */
+ double calcDistOnIntersecEdge(const WorldVector<double> &sP1,
+ const WorldVector<double> &sP2,
+ const WorldVector<double> &v,
+ double &lambda);
+
+ //calculates the distance between two points
+ double calc_dist_between_two_points(
+ const WorldVector<double> &point1,
+ const WorldVector<double> &point2);
+
+ void cases_for_lambda(const WorldVector<double> &SP,
+ const WorldVector<double> &SP_projected,
+ const WorldVector<double> &point1,
+ const WorldVector<double> &point2,
+ double &lambda,
+ double &dist);
+
+ //calculates the projection of the point "v" onto a straight line
+ //given by two points "sP1" nad "sP2"
+ void projected_on_a_straight_line(const WorldVector<double> &sP1,
+ const WorldVector<double> &sP2,
+ const WorldVector<double> &v,
+ WorldVector<double> &v_p,
+ double &lambda_out );
+
+ protected:
+ /**
+ * Object needed to extrapolate velocity from the interface.
+ */
+ VelocityExt *velExt;
+ };
+
+#endif // BOUNDARYELEMENTTOPDIST_H
diff --git a/AMDiS/Reinit/src/ElementLevelSet.cc b/AMDiS/Reinit/src/ElementLevelSet.cc
new file mode 100644
index 00000000..5ef6b076
--- /dev/null
+++ b/AMDiS/Reinit/src/ElementLevelSet.cc
@@ -0,0 +1,423 @@
+#include "ElementLevelSet.h"
+#include "ElInfo.h"
+
+int
+ElementLevelSet::createElementLevelSet(const ElInfo *elInfo_,
+ const bool doCalcIntersecPts_)
+{
+ Element *el = elInfo_->getElement();
+
+ if (elInfo == NULL || el != lastEl) {
+ /**
+ * Element has changed. New calculation.
+ */
+
+ // Set new element.
+ lastEl = el;
+
+ // Set information for element and reset data.
+ setElement(elInfo_);
+
+ // Calculate level set function values for each vertex of element.
+ calculateElementLevelSetVal();
+
+ // Check whether level set function values are not too small.
+ checkElementLevelSetVal();
+
+ // Calculate status of element and element vertices.
+ calculateElementStatus();
+
+ // Calculate intersection points with zero level set if necessary.
+ if (doCalcIntersecPts_ && elStatus == LEVEL_SET_BOUNDARY) {
+ calculateIntersecPoints();
+
+ // dim == 3: sort intersection points if necessary
+ if (dim == 3 && numIntersecPoints == 4)
+ sortIntersecPoints_4IP3D();
+ }
+ }
+
+// else {
+ /**
+ * LevelSet-Status for element has already been calculated.
+ */
+// }
+
+ return elStatus;
+}
+
+void
+ElementLevelSet::calculateElementLevelSetVal()
+{
+ DimVec<double> lambda(dim, DEFAULT_VALUE, 0.0);
+
+ for (int i = 0; i <= dim; i++) {
+
+ lambda[i] = 1.0;
+ lSFct->setElInfo(elInfo);
+ elVertexLevelSetVec[i] = (*lSFct)((const DimVec<double>) lambda);
+ lambda[i] = 0.0;
+ }
+}
+
+int
+ElementLevelSet::calculateElementStatus()
+{
+ for (int i=0; i<=dim; ++i) {
+ if (elVertexLevelSetVec[i] < 0) {
+ elVertexStatusVec[i] = LEVEL_SET_INTERIOR;
+ numElVertexInterior++;
+ }
+ else {
+ elVertexStatusVec[i] = LEVEL_SET_EXTERIOR;
+ numElVertexExterior++;
+ }
+ }
+
+ /**
+ * Calculate level set status of element.
+ */
+ if (numElVertexInterior == dim+1) {
+ elStatus = LEVEL_SET_INTERIOR;
+ return LEVEL_SET_INTERIOR;
+ }
+ if (numElVertexExterior == dim+1) {
+ elStatus = LEVEL_SET_EXTERIOR;
+ return LEVEL_SET_EXTERIOR;
+ }
+ elStatus = LEVEL_SET_BOUNDARY;
+ return LEVEL_SET_BOUNDARY;
+}
+
+void
+ElementLevelSet::calculateIntersecPoints()
+{
+ //**************************************************************************
+ // The level set function is linearly approximated by a hyperplane through
+ // the points of the graph of the level set function in the vertices
+ // of the element.
+ // This routine calculates the intersection points of the hyperplane
+ // with the edges of the element.
+ //**************************************************************************
+ DimVec<double> tempPoint(dim);
+ DimVec<double> zeroVec(dim, DEFAULT_VALUE, 0.0);
+ int i, j;
+
+ /**
+ * Get intersection points (in barycentric coordinates).
+ *
+ * An edge is intersected if its vertices have level sets with
+ * different sign.
+ */
+ for (i = 0; i <= dim; i++) {
+ for (j = i+1; j <= dim; j++) {
+
+ if (elVertexStatusVec[i] * elVertexStatusVec[j] < 0.0) {
+
+ tempPoint = zeroVec;
+ tempPoint[j] = elVertexLevelSetVec[i] /
+ (elVertexLevelSetVec[i] - elVertexLevelSetVec[j]);
+ checkIntersecBary(tempPoint[j]);
+ tempPoint[i] = 1 - tempPoint[j];
+
+ (*elIntersecPoints)[numIntersecPoints] = tempPoint;
+ ++numIntersecPoints;
+ }
+
+ } // for(j ...
+ } // for(i ...
+}
+
+int
+ElementLevelSet::checkElementLevelSetVal()
+{
+ int changed = 0;
+ double abs_grd_val = 0.0;
+ double abs_min_val = 0.0;
+ double comp = 0.0;
+
+ for (int i=0; i<dim; ++i) {
+ comp = elVertexLevelSetVec[i]-elVertexLevelSetVec[dim];
+ abs_grd_val += comp*comp;
+ }
+ abs_grd_val = sqrt(abs_grd_val);
+
+ abs_min_val = LS_VAL_TOL*abs_grd_val;
+ abs_min_val = (abs_min_val > LS_VAL_MIN) ? abs_min_val : LS_VAL_MIN;
+
+ for (int i=0; i<=dim; ++i) {
+ if (fabs(elVertexLevelSetVec[i]) < abs_min_val) {
+// elVertexLevelSetVec[i] = abs_min_val;
+ elVertexLevelSetVec[i] = (elVertexLevelSetVec[i] < 0) ?
+ -abs_min_val : abs_min_val;
+ ++changed;
+ }
+ }
+
+ return changed;
+}
+
+bool
+ElementLevelSet::checkIntersecBary(double &bary)
+{
+ if (bary < SP_BARY_TOL) {
+ bary = SP_BARY_TOL;
+ return true;
+ }
+ if (bary > 1-SP_BARY_TOL) {
+ bary = 1-SP_BARY_TOL;
+ return true;
+ }
+
+ return false;
+}
+
+void
+ElementLevelSet::sortIntersecPoints_4IP3D()
+{
+ FUNCNAME("sortIntersecPoints_4IP3D");
+
+ int indexFace1 = 0;
+ int indexFace2 = 0;
+ int indexOpVert = 0;
+ DimVec<double> tempPoint(dim);
+ int i,j;
+
+ /**
+ * Consider the 4 intersection points S0, S1, S2 and S3. If the components
+ * of S0 indexFace1 and indexFace2 are not zero, i.e. S0 lies in
+ * the element faces indexFace1 and indexFace2, the intersection
+ * point with zero components indexFace1 and indexFace2 is
+ * opposite to S0 in the intersection plane. Move this vertex to the end
+ * of the intersection point list elIntersecPoints.
+ */
+
+ // Get indexFace1 and indexFace2.
+ for (i = 0; i < numIntersecPoints; i++) {
+ if (fabs((*elIntersecPoints)[0][i]) > 1.e-15) {
+ indexFace1 = i;
+ break;
+ }
+ }
+ for (j = i+1; j < numIntersecPoints; j++) {
+ if (fabs((*elIntersecPoints)[0][j]) > 1.e-15) {
+ indexFace2 = j;
+ break;
+ }
+ }
+
+ // Get index of vertex opposite to S0.
+ for (i = 1; i < numIntersecPoints; i++) {
+ if (fabs((*elIntersecPoints)[i][indexFace1]) <= 1.e-15 &&
+ fabs((*elIntersecPoints)[i][indexFace2]) <= 1.e-15) {
+ indexOpVert = i;
+ break;
+ }
+ }
+
+ // Move vertex to the end of \ref elIntersecPoints.
+ if (indexOpVert != numIntersecPoints-1) {
+ tempPoint = (*elIntersecPoints)[indexOpVert];
+ (*elIntersecPoints)[indexOpVert] = (*elIntersecPoints)[numIntersecPoints-1];
+ (*elIntersecPoints)[numIntersecPoints-1] = tempPoint;
+ }
+
+ return;
+}
+
+void
+ElementLevelSet::calcIntersecNormal(WorldVector<double> &normalVec)
+{
+ FUNCNAME("ElementLevelSet::calcIntersecNormal");
+
+ switch(dim) {
+ case 2: calcIntersecNormal_2d(normalVec);
+ break;
+ case 3: calcIntersecNormal_3d(normalVec);
+ break;
+ default: ERROR_EXIT("illegal dimension !\n");
+ }
+};
+
+void
+ElementLevelSet::calcIntersecNormal_2d(WorldVector<double> &normalVec)
+{
+ FUNCNAME("ElementLevelSet::calcIntersecNormal_2d");
+
+ TEST_EXIT(numIntersecPoints == 2)("illegal number of intersection points !\n");
+
+ // ===== Get world coordinates of intersection points. =====
+ WorldVector<double> sP1;
+ WorldVector<double> sP2;
+ elInfo->coordToWorld((*elIntersecPoints)[0], sP1);
+ elInfo->coordToWorld((*elIntersecPoints)[1], sP2);
+
+ // ===== Calculate normal vector. =====
+ double norm2 = 0.0;
+ double val;
+ for (int i=0; i<dim; ++i){
+ val = sP1[i] - sP2[i];
+ norm2 += val*val;
+ normalVec[dim-1-i] = val;
+ }
+ normalVec[0] *= -1;
+ double norm = sqrt(norm2);
+ for(int i=0; i<dim; ++i) {
+ normalVec[i] = 1/norm * normalVec[i];
+ }
+
+ // ===== Set correct orientation (exterior normal vector). =====
+
+ // Calculate barycenter.
+ WorldVector<double> baryc;
+ for (int i=0; i<dim; ++i) {
+ baryc[i] = 0.5*sP1[i] + 0.5*sP2[i];
+ }
+ double d = sqrt((sP1-sP2)*(sP1-sP2));
+
+ // Barycenter + factor*normalVec.
+ double elSize = sqrt(fabs(elInfo->getDet()));
+ double factor = 0.01*d/elSize;
+ WorldVector<double> tmpPoint;
+ int cntr = 0;
+ DimVec<double> lambda(dim, NO_INIT);
+
+ while (1) {
+ ++cntr;
+
+ for(int i=0; i<dim; ++i) {
+ tmpPoint[i] = baryc[i] + factor*normalVec[i];
+ }
+
+ elInfo->worldToCoord(tmpPoint, &lambda);
+ for (int i=0; i<=dim; ++i) {
+ if (lambda[i] < 0) {
+ factor *= 0.1;
+
+ if (cntr == 10) {
+ WARNING("inefficient normal vector calculation !\n");
+ }
+ if (cntr == 100) {
+ ERROR_EXIT("infinite loop !\n");
+ }
+
+ continue;
+ }
+ }
+
+ break;
+ }
+
+ if (ElementLevelSet::calcLevelSetFct(lambda) < 0) {
+ for (int i=0; i<dim; ++i) {
+ normalVec[i] *= -1;
+ }
+ }
+};
+
+void
+ElementLevelSet::calcIntersecNormal_3d(WorldVector<double> &normalVec)
+{
+ FUNCNAME("ElementLevelSet::calcIntersecNormal_3d");
+
+ TEST_EXIT(numIntersecPoints == 3 || numIntersecPoints == 4)("illegal number of intersection points !\n");
+
+ // ===== Get world coordinates of intersection points. =====
+ WorldVector<double> sP1;
+ WorldVector<double> sP2;
+ WorldVector<double> sP3;
+ elInfo->coordToWorld((*elIntersecPoints)[0], sP1);
+ elInfo->coordToWorld((*elIntersecPoints)[1], sP2);
+ elInfo->coordToWorld((*elIntersecPoints)[2], sP3);
+
+ // ===== Calculate normal vector. =====
+ WorldVector<double> A = sP2-sP1;
+ WorldVector<double> B = sP3-sP1;
+ vectorProduct(A, B, normalVec);
+
+ double norm = sqrt(normalVec * normalVec);
+ for(int i=0; i<dim; ++i) {
+ normalVec[i] = 1/(norm+1.e-4) * normalVec[i];
+ }
+
+ // ===== Set correct orientation (exterior normal vector). =====
+
+ // Calculate barycenter.
+ WorldVector<double> baryc;
+ for (int i=0; i<dim; ++i) {
+ baryc[i] = 1.0/3*sP1[i] + 1.0/3*sP2[i] + 1.0/3*sP3[i];
+ }
+ double d = sqrt((sP1-sP2)*(sP1-sP2));
+
+ // Barycenter + factor*normalVec.
+ double elSize = pow(fabs(elInfo->getDet()), 1.0/3);
+ double factor = 0.01*d/elSize;
+ WorldVector<double> tmpPoint;
+ int cntr = 0;
+ DimVec<double> lambda(dim, NO_INIT);
+
+ while (1) {
+ ++cntr;
+
+ for(int i=0; i<dim; ++i) {
+ tmpPoint[i] = baryc[i] + factor*normalVec[i];
+ }
+
+ elInfo->worldToCoord(tmpPoint, &lambda);
+ for (int i=0; i<=dim; ++i) {
+ if (lambda[i] < 0) {
+ factor *= 0.1;
+
+ if (cntr == 10) {
+ WARNING("inefficient normal vector calculation !\n");
+ }
+ if (cntr == 100) {
+ ERROR_EXIT("infinite loop !\n");
+ }
+
+ continue;
+ }
+ }
+
+ break;
+ }
+
+ if (ElementLevelSet::calcLevelSetFct(lambda) < 0) {
+ for (int i=0; i<dim; ++i) {
+ normalVec[i] *= -1;
+ }
+ }
+};
+
+int
+ElementLevelSet::getVertexPos(const DimVec<double> barCoords)
+{
+ double vertex_val;
+
+ for (int i=0; i<=dim; ++i) {
+ if (barCoords[i] > 1-1.e-15) {
+ vertex_val = elVertexLevelSetVec[i];
+ break;
+ }
+ }
+
+ if (vertex_val > 0) {
+ return LEVEL_SET_EXTERIOR;
+ }
+ else {
+ return LEVEL_SET_INTERIOR;
+ }
+}
+
+int
+ElementLevelSet::getElPos(const DimVec<double> barCoords)
+{
+ double ls_val = calcLevelSetFct(barCoords);
+ if (ls_val > 0) {
+ return LEVEL_SET_EXTERIOR;
+ }
+ else {
+ return LEVEL_SET_INTERIOR;
+ }
+}
+
diff --git a/AMDiS/Reinit/src/ElementLevelSet.h b/AMDiS/Reinit/src/ElementLevelSet.h
new file mode 100644
index 00000000..5bcf7964
--- /dev/null
+++ b/AMDiS/Reinit/src/ElementLevelSet.h
@@ -0,0 +1,497 @@
+#ifndef AMDIS_ELEMENTLEVELSET_H
+#define AMDIS_ELEMENTLEVELSET_H
+
+#include "ElementFunction.h"
+#include "FixVec.h"
+// #include "MemoryManager.h"
+#include "Parameters.h"
+
+namespace AMDiS {
+ class Element;
+ class ElInfo;
+ class Mesh;
+}
+
+using namespace AMDiS;
+using namespace std;
+
+// ===========================================================================
+// ===== class ElementLevelSet ===============================================
+// ===========================================================================
+//
+// Class Description:
+// The class ElementLevelSet contains the functionality for
+// - calculating the intersection points resulting from an intersection
+// of the boundary (level set zero) with an element,
+// - calculating the status of an element, i.e. is the element cut
+// by the zero level set or not.
+// The calculation of the intersection points is done as follows:
+// The level set function is linearly approximated on the element, i.e. its
+// graph is approximated by the plane through the level set values of the
+// element vertices. We approximate the intersection points by the
+// intersection of the plane with the element edges.
+// If in 3-dimensional finite element space the intersection produced 4
+// intersection points S0, S1, S2, S3, the intersection points in
+// elIntersecPoints are rearranged so that S1 and S2 divides the intersection
+// plane in two (dim - 1)-dimensional simplices.
+//
+// Constants indicating the level set status of element:
+// LEVEL_SET_INTERIOR - element is in domain where level set function
+// is negative
+// LEVEL_SET_BOUNDARY - element is in domain where level set function
+// is positive
+// LEVEL_SET_EXTERIOR - element is cut by the zero level set
+//
+// Main routines:
+// setElementLevelSet() - Defines the level set function for the
+// following calculations.
+// createElementLevelSet() - Calculates level set status of element and
+// intersection points if needed.
+// calculateElementLevelSetVal() - Calculates values of the level set
+// function in the element vertices.
+// setElement() - Sets elInfo.
+// getElementLevelSetStatus() - Returns status of element.
+// getElementIntersecPoints() - Returns intersection points.
+// getElVertStatusVec() - Returns vector with status information
+// for each vertex.
+// ===========================================================================
+class ElementLevelSet
+{
+ public:
+ MEMORY_MANAGED(ElementLevelSet);
+
+ /**
+ * Constructor.
+ */
+ ElementLevelSet(const char *name_,
+ ElementFunction<double> *lSFct_,
+ Mesh *mesh_)
+ : name(name_),
+ elInfo(NULL),
+ lastEl(NULL),
+ level_set_domain(LEVEL_SET_UNDEFINED),
+ numIntersecPoints(0),
+ elStatus(LEVEL_SET_UNDEFINED),
+ numElVertexInterior(0),
+ numElVertexBoundary(0),
+ numElVertexExterior(0),
+ LS_VAL_TOL(1.e-8),
+ LS_VAL_MIN(1.e-8),
+ SP_BARY_TOL(1.e-7)
+ {
+ FUNCNAME("ElementLevelSet::ElementLevelSet()");
+
+ TEST_EXIT(lSFct_ || mesh_)
+ ("illegal initialization of ElementLevelSet!\n");
+
+ lSFct = lSFct_;
+ mesh = mesh_;
+ dim = mesh->getDim();
+
+ elIntersecPoints =
+ NEW VectorOfFixVecs<DimVec<double> >(dim,
+ MAX_INTERSECTION_POINTS,
+ NO_INIT);
+ elVertexStatusVec = new int[dim+1];
+ elVertexLevelSetVec = new double[dim+1];
+
+ int setElementLevelSetTol = 0;
+ GET_PARAMETER(0, name + "->set ElementLevelSet tolerances", "%d",
+ &setElementLevelSetTol);
+ if (setElementLevelSetTol) {
+
+ GET_PARAMETER(0, name + "->LS_VAL_TOL", "%f", &LS_VAL_TOL);
+ GET_PARAMETER(0, name + "->LS_VAL_MIN", "%f", &LS_VAL_MIN);
+ GET_PARAMETER(0, name + "->SP_BARY_TOL", "%f", &SP_BARY_TOL);
+
+ TEST_EXIT(LS_VAL_TOL > 0)("illegal LS_VAL_TOL\n");
+ TEST_EXIT(LS_VAL_MIN > 0)("illegal LS_VAL_MIN\n");
+ TEST_EXIT(SP_BARY_TOL > 0)("illegal SP_BARY_TOL\n");
+ }
+ };
+
+ /**
+ * Destructor.
+ */
+ ~ElementLevelSet()
+ {
+ if (elVertexStatusVec)
+ delete [] elVertexStatusVec;
+ if(elVertexLevelSetVec)
+ delete [] elVertexLevelSetVec;
+ if (elIntersecPoints)
+ DELETE elIntersecPoints;
+ };
+
+ /**
+ * Calculates LevelSet-status of element and its intersection points
+ * with the zero level set if necessary.
+ *
+ * Result:
+ * LEVEL_SET_BOUNDARY: Element elInfo is intersected by levelSetFct.
+ * LEVEL_SET_EXTERIOR / LEVEL_SET_INTERIOR: Element lies completely on
+ * one side of the zero level set.
+ *
+ * After proceeding this function, information about the level set
+ * status is given in:
+ * elStatus: status of element (LEVEL_SET_BOUNDARY, LEVEL_SET_INTERIOR or
+ * EXTERIOR)
+ * elVertexStatusVec: stores status of each vertex of element
+ * elVertexLevelSetVec: stores level set function value of each vertex of
+ * element
+ * numElVertexInterior: number of vertices of element with status
+ * LEVEL_SET_INTERIOR
+ * numElVertexExterior: number of vertices of element with status
+ * LEVEL_SET_EXTERIOR
+ * numElVertexBoundary: number of vertices of element with status
+ * LEVEL_SET_BOUNDARY
+ * elIntersecPoints: stores the intersection points produced by the
+ * intersection of element with the zero level set
+ * numIntersecPoints: number of intersection points
+ */
+ int createElementLevelSet(const ElInfo *elInfo_,
+ const bool doCalcIntersecPts_ = true);
+
+ /**
+ * Gets value of level set function at point given in
+ * barycentric coordinates.
+ */
+ inline double calcLevelSetFct(const DimVec<double>& bary) {
+ return (*lSFct)(bary);
+ };
+
+ /**
+ * Resets level set information on element.
+ */
+ inline void resetElement() {
+ FUNCNAME("ElementLevelSet::resetElement");
+
+ numElVertexInterior = 0;
+ numElVertexBoundary = 0;
+ numElVertexExterior = 0;
+ numIntersecPoints = 0;
+ elStatus = LEVEL_SET_UNDEFINED;
+ };
+
+ /**
+ * Defines current element (elInfo).
+ */
+ inline void setElement(const ElInfo *elInfo_) {
+ elInfo = elInfo_;
+ resetElement();
+ };
+
+ /**
+ * Set level_set_domain.
+ */
+ inline void setLevelSetDomain(int status_) {
+
+ TEST_EXIT(status_ == LEVEL_SET_INTERIOR ||
+ status_ == LEVEL_SET_EXTERIOR ||
+ status_ == LEVEL_SET_BOUNDARY)("illegal level set status !\n");
+ level_set_domain = status_;
+ };
+
+ /**
+ * Functions to set tolerances for intersection point calculation.
+ */
+ inline void setLsValTol(double tol) {LS_VAL_TOL = tol;};
+ inline void setLsValMin(double min) {LS_VAL_MIN = min;};
+ inline void setSpBaryTol(double tol) {SP_BARY_TOL = tol;};
+
+ /**
+ * Get level_set_domain.
+ */
+ inline const int& getLevelSetDomain() const {
+ return level_set_domain;
+ };
+
+ /**
+ * Get LevelSet-Status of element.
+ */
+ inline const int& getElementLevelSetStatus() const {
+ return elStatus;
+ };
+
+ /**
+ * Get number of vertices which are intersection points.
+ */
+ inline const int& getNumVertIntPoints() const {
+ FUNCNAME("ElementLevelSet::getNumVertIntPoints");
+ TEST_EXIT(numElVertexBoundary == 0)("numElVertexBoundary should be zero!\n");
+ return numElVertexBoundary;
+ };
+
+ /**
+ * Get vector elVertexStatusVec.
+ */
+ inline const int *getElVertStatusVec() const {
+ return elVertexStatusVec;
+ };
+
+ /**
+ * Get i-th component of vector elVertexLevelSetVec.
+ */
+ inline const double getElVertLevelSetVec(const int i) const {
+ return elVertexLevelSetVec[i];
+ };
+
+ /**
+ * Get vector elVertexLevelSetVec.
+ */
+ inline const double *getElVertLevelSetVec() const {
+ return elVertexLevelSetVec;
+ };
+
+ /**
+ * Get levelSetFct.
+ */
+ inline ElementFunction<double> *getLevelSetFct() const {
+ return lSFct;
+ };
+
+ /**
+ * Get mesh.
+ */
+ inline Mesh *getMesh() const {
+ return mesh;
+ };
+
+ /**
+ * Get dim.
+ */
+ inline int getDim() const {
+ return dim;
+ };
+
+ /**
+ * Get the intersection points.
+ */
+ inline VectorOfFixVecs<DimVec<double> > *getElIntersecPoints() const {
+ return elIntersecPoints;
+ };
+
+ /**
+ * Get number of intersection points.
+ */
+ inline int getNumElIntersecPoints() const {
+ return numIntersecPoints;
+ };
+
+ /**
+ * Calculate exterior normal to intersection plane.
+ */
+ void calcIntersecNormal(WorldVector<double> &normal);
+
+ /**
+ * Gets position of point in element with barycentric coordinates
+ * barCoords, i.e. whether point is in the domain with positive
+ * (LEVEL_SET_EXTERIOR) or negative (LEVEL_SET_INTERIOR) level set
+ * function values. Uses level set function, thus element vertices
+ * may have level set function value zero.
+ */
+ int getElPos(const DimVec<double> barCoords);
+
+ /**
+ * Gets position of element vertex given in barycentric coordinates,
+ * i.e. whether element vertex is in the domain with positive
+ * (LEVEL_SET_EXTERIOR) or negative (LEVEL_SET_INTERIOR) level set
+ * function values. Uses elVertexLevelSetVec.
+ */
+ int getVertexPos(const DimVec<double> barCoords);
+
+ protected:
+ /**
+ * Calculates level set value of each vertex of element.
+ */
+ void calculateElementLevelSetVal();
+
+ /**
+ * Calculates the status of an element.
+ *
+ * Note: Uses values in elVertexLevelSetVec.
+ *
+ * Return value:
+ * LEVEL_SET_INTERIOR element lies completely inside
+ * LEVEL_SET_EXTERIOR element lies completely outside
+ * LEVEL_SET_BOUNDARY element is cut by the zero level set
+ */
+ int calculateElementStatus();
+
+ /**
+ * Calculates intersection points of zero level set with element.
+ *
+ * Note: Uses elVertexLevelSet.
+ */
+ void calculateIntersecPoints();
+
+ /**
+ * Checks whether level set values of element (in elVertexLevelSetVec)
+ * are below a certain bound and corrects them if this is the case.
+ *
+ * Return value: number of values corrected
+ */
+ int checkElementLevelSetVal();
+
+ /**
+ * Checks whether barycentric coordinate of intersection point is not
+ * too small and corrects it if this is the case.
+ *
+ * Return value: true - barycentric coordinate has been corrected
+ * false - barycentric coordinate is ok
+ */
+ bool checkIntersecBary(double &bary);
+
+ /**
+ * Sort intersection points S0, S1, S2 and S3 in \ref elIntersecPoints in
+ * such a way that afterwards, a line through S1 and S2 divides the
+ * intersection plane into two (\ref dim - 1)-dimensional simplices.
+ */
+ void sortIntersecPoints_4IP3D();
+
+ /**
+ * Calculate exterior normal to intersection plane for dimension 2.
+ */
+ void calcIntersecNormal_2d(WorldVector<double> &normal);
+
+ /**
+ * Calculate exterior normal to intersection plane for dimension 3.
+ */
+ void calcIntersecNormal_3d(WorldVector<double> &normal);
+
+ public:
+ /**
+ * Constants characterizing element position with respect to zero level set.
+ */
+ static const int LEVEL_SET_INTERIOR = -1;
+ static const int LEVEL_SET_BOUNDARY = 0;
+ static const int LEVEL_SET_EXTERIOR = 1;
+ static const int LEVEL_SET_UNDEFINED = -2;
+
+ protected:
+ /**
+ * Name of this object.
+ */
+ std::string name;
+
+ /**
+ * Level set function.
+ */
+ ElementFunction<double> *lSFct;
+
+ /**
+ * Mesh.
+ */
+ Mesh *mesh;
+
+ /**
+ * elInfo of element.
+ */
+ const ElInfo *elInfo;
+
+ /**
+ * Pointer to last element processed calculations on whithin this class.
+ */
+ Element *lastEl;
+
+ /**
+ * Indicator which can be used for example for function evaluation
+ * or integration on subelements. Indicates whether point/subelement ...
+ * is inside (LEVEL_SET_INTERIOR) or outside (LEVEL_SET_EXTERIOR) the
+ * zero level set or is cut by the zero level set (LEVEL_SET_BOUNDARY).
+ */
+ int level_set_domain;
+
+ /**
+ * Dimension of the problem. dim + 1 is the number of vertices
+ * of element.
+ */
+ int dim;
+
+ /**
+ * Vector for intersection points produced by the intersection of linearly
+ * approximated level set function with the edges of element.
+ */
+ VectorOfFixVecs<DimVec<double> > *elIntersecPoints;
+
+ /**
+ * Number of intersection points.
+ */
+ int numIntersecPoints;
+
+ /**
+ * LevelSet-Status of element.
+ */
+ int elStatus;
+
+ /**
+ * Holds for each vertex of element the information about the position
+ * of the vertex with respect to the zero level set.
+ */
+ int *elVertexStatusVec;
+
+ /**
+ * Stores for each vertex of element the level set of the vertex.
+ */
+ double *elVertexLevelSetVec;
+
+ /**
+ * Number of vertices in element with level set status LEVEL_SET_INTERIOR.
+ */
+ int numElVertexInterior;
+
+ /**
+ * Number of vertices in element with level set status LEVEL_SET_BOUNDARY.
+ *
+ * Note: should be zero
+ */
+ int numElVertexBoundary;
+
+ /**
+ * Number of vertices in element with level set status LEVEL_SET_EXTERIOR.
+ */
+ int numElVertexExterior;
+
+ /**
+ * Tolerance used in the calculation of the local minimal level set value.
+ * The local minimal level set value depends on the gradient of the
+ * level set function.
+ * Used for the calculation of intersection points.
+ *
+ * If intersection points are too close to vertices, they are slightly
+ * moved.
+ * IDEA: If d is the distance of an intersection point to vertex v,
+ * the property
+ *
+ * d > LS_VAL_TOL * h
+ *
+ * must be satisfied.
+ * In the implementation this results in
+ *
+ * phi(v) > LS_VAL_TOL * h * grad .
+ */
+ double LS_VAL_TOL;
+
+ /**
+ * Lower bound for level set value on elements cut by the zero level set.
+ * Used for the calculation of intersection points.
+ */
+ double LS_VAL_MIN;
+
+ /**
+ * Lower bound for barycentric coordinates of intersection points.
+ *
+ * Each component x of the barycentric coordinates of an intersection
+ * point satisfies
+ *
+ * SP_BARY_TOL < x < 1 - SP_BARY_TOL .
+ */
+ double SP_BARY_TOL;
+
+ /*
+ * Maximum number of intersection points.
+ */
+ static const int MAX_INTERSECTION_POINTS = 4;
+};
+
+#endif // AMDIS_ELEMENTLEVELSET_H
diff --git a/AMDiS/Reinit/src/ElementUpdate.h b/AMDiS/Reinit/src/ElementUpdate.h
new file mode 100644
index 00000000..c20fb738
--- /dev/null
+++ b/AMDiS/Reinit/src/ElementUpdate.h
@@ -0,0 +1,43 @@
+#ifndef ELEMENTUPDATE_H
+#define ELEMENTUPDATE_H
+
+#include "FixVec.h"
+// #include "MemoryManager.h"
+
+#include "VelocityExt.h"
+
+using namespace AMDiS;
+
+class ElementUpdate
+{
+ public:
+ MEMORY_MANAGED(ElementUpdate);
+
+ /**
+ * Constructor.
+ */
+ ElementUpdate(VelocityExt *velExt_)
+ : velExt(velExt_)
+ {};
+
+ /**
+ * Virtual destructor.
+ */
+ virtual ~ElementUpdate() {};
+
+ /**
+ * Pure virtual function.
+ * Calculate Bornemann update on an element. Realization is implemented for
+ * dimensions 2 and 3.
+ */
+ virtual double calcElementUpdate(const FixVec<WorldVector<double> *, VERTEX> &vert,
+ FixVec<double, VERTEX> &uhVal) = 0;
+
+ protected:
+ /**
+ * Object needed to extrapolate velocity from the interface.
+ */
+ VelocityExt *velExt;
+};
+
+#endif // ELEMENTUPDATE_H
diff --git a/AMDiS/Reinit/src/ElementUpdate_2d.cc b/AMDiS/Reinit/src/ElementUpdate_2d.cc
new file mode 100644
index 00000000..92d76aee
--- /dev/null
+++ b/AMDiS/Reinit/src/ElementUpdate_2d.cc
@@ -0,0 +1,52 @@
+#include "ElementUpdate_2d.h"
+
+double
+ElementUpdate_2d::calcElementUpdate(
+ const FixVec<WorldVector<double> *, VERTEX> &vert,
+ FixVec<double, VERTEX> &uhVal)
+{
+ WorldVector<double> xhminusYh = *(vert[2]) - *(vert[0]);
+ WorldVector<double> zhminusYh = *(vert[1]) - *(vert[0]);
+ WorldVector<double> xhminusZh = *(vert[2]) - *(vert[1]);
+
+ double norm_zhminusYh = sqrt(zhminusYh * zhminusYh);
+ double norm_xhminusYh = sqrt(xhminusYh * xhminusYh);
+ double norm_xhminusZh = sqrt(xhminusZh * xhminusZh);
+
+ double delta = (uhVal[1]-uhVal[0]) / norm_zhminusYh;
+
+ double c_alpha = 0;
+ double c_beta = 0;
+ double sP = xhminusYh * zhminusYh;
+ c_alpha = sP / (norm_xhminusYh * norm_zhminusYh);
+ sP = xhminusZh * zhminusYh;
+ c_beta = -sP / (norm_xhminusZh * norm_zhminusYh);
+
+ double update = 0;
+ if (c_alpha <= delta) {
+ update = uhVal[0] + norm_xhminusYh;
+ //save barycentric coordinates for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_2D(1,0,0);
+ }
+ } else if (delta <= -c_beta) {
+ update = uhVal[1] + norm_xhminusZh;
+ //save barycentric coordinates for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_2D(0,1,0);
+ }
+ } else {
+ update = uhVal[0] +
+ (c_alpha*delta + sqrt((1-c_alpha*c_alpha)*(1-delta*delta))) *
+ norm_xhminusYh;
+ //calculate and save barycentric coordinates for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->calcBarycentricCoords_2D(delta, c_alpha, norm_zhminusYh, norm_xhminusYh);
+ }
+ }
+
+ return update;
+}
diff --git a/AMDiS/Reinit/src/ElementUpdate_2d.h b/AMDiS/Reinit/src/ElementUpdate_2d.h
new file mode 100644
index 00000000..a5fbbd96
--- /dev/null
+++ b/AMDiS/Reinit/src/ElementUpdate_2d.h
@@ -0,0 +1,33 @@
+#ifndef ELEMENTUPDATE_2D_H
+#define ELEMENTUPDATE_2D_H
+
+#include "FixVec.h"
+// #include "MemoryManager.h"
+
+#include "ElementUpdate.h"
+#include "VelocityExt.h"
+
+using namespace AMDiS;
+
+class ElementUpdate_2d : public ElementUpdate
+{
+ public:
+ MEMORY_MANAGED(ElementUpdate_2d);
+
+ /**
+ * Constructor.
+ */
+ ElementUpdate_2d(VelocityExt *velExt_ = NULL)
+ : ElementUpdate(velExt_)
+ {};
+
+ /**
+ * Realization of ElementUpdate::calcElementUpdate.
+ * Calculates the Bornemann update for an element.
+ */
+ double calcElementUpdate(const FixVec<WorldVector<double> *, VERTEX> &vert,
+ FixVec<double, VERTEX> &uhVal);
+};
+
+#endif // ELEMENTUPDATE_2D_H
+
diff --git a/AMDiS/Reinit/src/ElementUpdate_3d.cc b/AMDiS/Reinit/src/ElementUpdate_3d.cc
new file mode 100644
index 00000000..f21bef58
--- /dev/null
+++ b/AMDiS/Reinit/src/ElementUpdate_3d.cc
@@ -0,0 +1,609 @@
+#include "ElementUpdate_3d.h"
+
+double
+ElementUpdate_3d::calcElementUpdate(
+ const FixVec<WorldVector<double> *, VERTEX> &vert_,
+ FixVec<double, VERTEX> &uhVal)
+{
+ FUNCNAME("ElementUpdate_3d::calcElementUpdate");
+
+ int dim = 3;
+ double update = 0.0;
+ WorldVector<double> B1;
+ WorldVector<double> C1;
+ WorldVector<double> D1;
+ double uhVal0_orig; // original value of vertex 0 after sort and
+ // translation of vertices but before
+ // translation of function values to zero
+ FixVec<WorldVector<double> *, VERTEX> vert(dim, NO_INIT);
+ for (int i=0; i<=dim; ++i) {
+ vert[i] = vert_[i];
+ }
+
+ // ===== Sort vertices and translate element to point of origin. =====
+ sortAndTranslateElement(vert, uhVal, uhVal0_orig);
+
+ // ===== Rotate element.
+ // The rotated element is stored in B1, C1, D1 (A = zero vector
+ // after sort and translation and is fix under rotation).
+ // Return value: true - rotation exists
+ // false - rotation does not exist
+ bool rotationFlag = rotateElement(vert, uhVal, B1, C1, D1);
+
+ // ===== rotation exists ==> calculate projection and update value
+ // depending on the position of the
+ // projection point
+ // rotation does not exist ==> calculate update values on the
+ // faces of element
+
+ if (rotationFlag) {
+
+ // Calculate projection of B1, C1 and D1 on x-y-plane -> B2, C2, D2.
+ WorldVector<double> B2;
+ WorldVector<double> C2;
+ WorldVector<double> D2;
+
+ B2[0] = B1[0]; B2[1] = B1[1]; B2[2] = 0.0;
+ C2[0] = C1[0]; C2[1] = C1[1]; C2[2] = 0.0;
+ D2[0] = D1[0]; D2[1] = D1[1]; D2[2] = 0.0;
+
+ // Calculate position of D2 with respect to A, B2 and C2.
+ // If D2 lies in the triangle AB2C2, this is the minimizing point
+ // of the update formula. Else, the minimizing point lies on edges
+ // of this triangle and the 2-dimensional update formula is applied
+ // on faces of element to get the minimizing point.
+ int posFlag = calcPosition(B2, C2, D2);
+
+ WorldVector<double> tmpVec;
+ switch(posFlag) {
+
+ case IN_ABC: // update value is the distance of D1 to the x-y-plane
+ update = D1[2];
+ break;
+
+ case VERT_A: // minimizing point is vertex A
+ update = uhVal[0] + sqrt(*(vert[3]) * *(vert[3]));
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_3D(1,0,0,0);
+ }
+ break;
+
+ case VERT_B: // minimizing point is vertex B
+ tmpVec = (*(vert[3])) - (*(vert[1]));
+ update = uhVal[1] + sqrt(tmpVec * tmpVec);
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_3D(0,1,0,0);
+ }
+ break;
+
+ case VERT_C: // minimizing point is vertex C
+ tmpVec = (*(vert[3])) - (*(vert[2]));
+ update = uhVal[2] + sqrt(tmpVec * tmpVec);
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_3D(0,0,1,0);
+ }
+ break;
+
+ case EDGE_AB: // minimizing point lies on edge AB
+ update = calcFaceUpdate(vert[0], vert[1], vert[3],
+ uhVal[0], uhVal[1]);
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->copyAndExpandFaceCoords_3D(2,0);
+ }
+ break;
+
+ case EDGE_AC: // minimizing point lies on edge AC
+ update = calcFaceUpdate(vert[0], vert[2], vert[3],
+ uhVal[0], uhVal[2]);
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->copyAndExpandFaceCoords_3D(1,0);
+ }
+ break;
+
+ case EDGE_BC: // minimizing point lies on edge BC
+ update = calcFaceUpdate(vert[1], vert[2], vert[3],
+ uhVal[1], uhVal[2]);
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->copyAndExpandFaceCoords_3D(0,0);
+ }
+ break;
+
+ default: ERROR_EXIT("illegal position flag !\n");
+ break;
+ }
+
+ }
+ else { // rotation does not exist
+
+ // ===== Calculate the 2-dimensional update for D on all faces
+ // adjacent with D. The update value is the minimum of these
+ // updates. =====
+ double tmp_update;
+
+ update = calcFaceUpdate(vert[0], vert[1], vert[3],
+ uhVal[0], uhVal[1]);
+ //set barycentric coordinats for calculation of the velocity
+ //save index of element face
+ if(velExt != NULL)
+ {
+ velExt->copyAndExpandFaceCoords_3D(2,1);
+ velExt->setIndexFaceWithShortestDist(1);
+ }
+ tmp_update = calcFaceUpdate(vert[0], vert[2], vert[3],
+ uhVal[0], uhVal[2]);
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->copyAndExpandFaceCoords_3D(1,2);
+ }
+ if (tmp_update < update) {
+ update = tmp_update;
+ //save index of element face, if diatance is shorter
+ if(velExt != NULL)
+ {
+ velExt->setIndexFaceWithShortestDist(2);
+ }
+ }
+ tmp_update = calcFaceUpdate(vert[1], vert[2], vert[3],
+ uhVal[1], uhVal[2]);
+ //set barycentric coordinats for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->copyAndExpandFaceCoords_3D(0,3);
+ }
+ if (tmp_update < update) {
+ update = tmp_update;
+ //save index of element face, if distance is shorter
+ if(velExt != NULL)
+ {
+ velExt->setIndexFaceWithShortestDist(3);
+ }
+ }
+ }
+
+ // ===== Correct initial translation of function values to zero. =====
+ update += uhVal0_orig;
+
+ return update;
+}
+
+void
+ElementUpdate_3d::sortAndTranslateElement(
+ FixVec<WorldVector<double> *, VERTEX> &vert,
+ FixVec<double, VERTEX> &uhVal,
+ double &uhVal0_orig)
+{
+ WorldVector<double> *tmpVec;
+ double tmpVal;
+
+ // ===== Sort vertices. =====
+ if (uhVal[0] > uhVal[1]) {
+
+ if (uhVal[1] > uhVal[2]) {
+
+ tmpVec = vert[0];
+ vert[0] = vert[2];
+ vert[2] = tmpVec;
+
+ tmpVal = uhVal[0];
+ uhVal[0] = uhVal[2];
+ uhVal[2] = tmpVal;
+ //save permutation of vertices
+ if(velExt != NULL)
+ {
+ velExt->swapVertices(0,2);
+ }
+ }
+ else if (uhVal[0] > uhVal[2]) {
+
+ tmpVec = vert[0];
+ vert[0] = vert[1];
+ vert[1] = vert[2];
+ vert[2] = tmpVec;
+
+ tmpVal = uhVal[0];
+ uhVal[0] = uhVal[1];
+ uhVal[1] = uhVal[2];
+ uhVal[2] = tmpVal;
+ //save permutation of vertices
+ if(velExt != NULL)
+ {
+ velExt->swapVertices(0,1);
+ velExt->swapVertices(2,1);
+ }
+ }
+ else {
+
+ tmpVec = vert[0];
+ vert[0] = vert[1];
+ vert[1] = tmpVec;
+
+ tmpVal = uhVal[0];
+ uhVal[0] = uhVal[1];
+ uhVal[1] = tmpVal;
+ //save permutation of vertices
+ if(velExt != NULL)
+ {
+ velExt->swapVertices(0,1);
+ }
+ }
+ }
+ else if (uhVal[1] > uhVal[2]) {
+
+ if (uhVal[0] > uhVal[2]) {
+
+ tmpVec = vert[0];
+ vert[0] = vert[2];
+ vert[2] = vert[1];
+ vert[1] = tmpVec;
+
+ tmpVal = uhVal[0];
+ uhVal[0] = uhVal[2];
+ uhVal[2] = uhVal[1];
+ uhVal[1] = tmpVal;
+ //save permutation of vertices
+ if(velExt != NULL)
+ {
+ velExt->swapVertices(0,2);
+ velExt->swapVertices(1,2);
+ }
+ }
+ else {
+
+ tmpVec = vert[1];
+ vert[1] = vert[2];
+ vert[2] = tmpVec;
+
+ tmpVal = uhVal[1];
+ uhVal[1] = uhVal[2];
+ uhVal[2] = tmpVal;
+ //save permutation of vertices
+ if(velExt != NULL)
+ {
+ velExt->swapVertices(1,2);
+ }
+ }
+ }
+
+ // ===== Translation: A is moved to point of origin. =====
+ *(vert[1]) -= *(vert[0]);
+ *(vert[2]) -= *(vert[0]);
+ *(vert[3]) -= *(vert[0]);
+ *(vert[0]) = 0.0;
+
+ // ===== Adapt \ref uhVal. =====
+ uhVal[1] -= uhVal[0];
+ uhVal[2] -= uhVal[0];
+ uhVal0_orig = uhVal[0];
+ uhVal[0] = 0.0;
+}
+
+bool
+ElementUpdate_3d::rotateElement(
+ const FixVec<WorldVector<double> *, VERTEX> &vert,
+ const FixVec<double, VERTEX> &uhVal,
+ WorldVector<double> &B1,
+ WorldVector<double> &C1,
+ WorldVector<double> &D1)
+{
+ int orient = 1; // orientation of AB, AC and AD
+ WorldVector<double> B = *(vert[1]);
+ WorldVector<double> C = *(vert[2]);
+ WorldVector<double> D = *(vert[3]);
+ WorldVector<double> N; // unit normal for ABC
+ // (orientation of AB, AC and N is 1)
+ WorldVector<double> P; // projection of D on ABC
+ WorldVector<double> unit_B; // normalized AB
+ double norm_B; // norm of vector AB
+ double norm_C; // norm of vector AC
+ double norm_P; // norm of vector AP
+ double dD; // (signed) distance of D to ABC
+ double lambda, nue; // scalars
+ WorldVector<double> N1; // unit normal for A1B1C1
+ WorldVector<double> unit_B1; // normalized B1
+ double dB1 = uhVal[1]; // distance of B1 to x-y-plane
+ double dC1 = uhVal[2]; // distance of C1 to x-y-plane
+ WorldVector<double> E_21; // vector of rotated coordinate system
+ int dim = 3;
+ WorldVector<double> tmpVec;
+ double tmp;
+
+
+ // ===== Calculate orientation of ABCD. =====
+
+ // Get unit normal for ABC.
+ tmpVec[0] = B[1]*C[2] - B[2]*C[1];
+ tmpVec[1] = B[2]*C[0] - B[0]*C[2];
+ tmpVec[2] = B[0]*C[1] - B[1]*C[0];
+ tmp = sqrt(tmpVec * tmpVec);
+ for (int i=0; i<dim; ++i) { N[i] = tmpVec[i] / tmp; }
+
+ // Calculate distance of D to ABC.
+ dD = D * N;
+
+ // Orientation.
+ orient = (dD < 0) ? -1 : 1;
+
+
+ // ===== Calculate B1.
+ // We construct B1 in x-z-plane with distance dB1 to
+ // x-y-plane. =====
+
+ norm_B = sqrt(B * B);
+
+ tmp = norm_B*norm_B - dB1*dB1;
+ if (tmp < 1.e-15) {
+ return false; // rotation is not possible
+ }
+ B1[0] = sqrt(tmp);
+ B1[1] = 0.0;
+ B1[2] = dB1;
+
+
+ // ===== Calculate C1.
+ // If (E_1, E_2, E_3) is a coordinate system with
+ // AC = lambda*E_1 + nue*E_2, then for the rotated coordinate
+ // system (E_11, E_21, E_31) we have
+ // A1C1 = lambda*E_11 + nue*E_21.
+ //
+ // The coordinate system (E_1, E_2, E_3) does not need to
+ // be calculated, but we assume
+ // --> E_1 = unit_B,
+ // --> nue > 0. =====
+
+ // Calculate lambda and nue.
+ for (int i=0; i<dim; ++i) { unit_B[i] = B[i] / norm_B; }
+ norm_C = sqrt(C * C);
+
+ lambda = C * unit_B; // corresponds to norm(AC)*cos(AC,AB)
+ nue = sqrt(norm_C*norm_C - lambda*lambda);
+ // corresponds to norm(AC)*sin(AC,AB)
+
+ // Calculate E_11 (= unit_B1).
+ for (int i=0; i<dim; ++i) { unit_B1[i] = B1[i] / norm_B; }
+
+ // Calculate E_21.
+ E_21[2] = (dC1 - lambda*unit_B1[2])/nue;
+ // A1C1 = lambda*E_11 + nue*E_21 and C1[2] = dC1
+
+ E_21[0] = -unit_B1[2]*E_21[2]/unit_B1[0];
+ // E_11*E_21 = 0 and unit_B[1]=0
+
+ tmp = 1 - E_21[0]*E_21[0] - E_21[2]*E_21[2];
+ if (tmp < 1.e-15) {
+ return false; // rotation is not possible
+ }
+ E_21[1] = orient * sqrt(tmp); // norm(E_21) = 1 and
+ // in further construction: D1
+ // above x-y-plane
+
+ // Calculate C1.
+ C1[0] = lambda*unit_B1[0] + nue*E_21[0];
+ C1[1] = lambda*unit_B1[1] + nue*E_21[1];
+ C1[2] = lambda*unit_B1[2] + nue*E_21[2];
+
+
+ // ===== Calculate D1.
+ // If we have AP = lambda*E_1 + nue*E_2 for the projection P of
+ // D on ABC , then D1 = lambda*E_11 + nue*E_21 + dD*N1. =====
+
+ // Calculate P.
+ P[0] = D[0] - dD*N[0];
+ P[1] = D[1] - dD*N[1];
+ P[2] = D[2] - dD*N[2];
+
+ // Calculate lambda and nue.
+ norm_P = sqrt(P * P);
+ lambda = P * unit_B; // corresponds to norm(AP)*cos(AP,AB)
+ nue = sqrt(norm_P*norm_P - lambda*lambda);
+ // corresponds to norm(AP)*sin(AP,AB)
+ // Fix problem for angle(E_2,P) > Pi/2:
+ // project P onto E_2-axes (-> tmpVec),
+ // if angle(tmpVec,AC) > Pi/2 the scalar nue has negative sign.
+ tmpVec = P - unit_B*lambda;
+ if (tmpVec * C < 0) { nue *= -1; }
+
+ // Calculate N1.
+ tmpVec[0] = B1[1]*C1[2] - B1[2]*C1[1];
+ tmpVec[1] = B1[2]*C1[0] - B1[0]*C1[2];
+ tmpVec[2] = B1[0]*C1[1] - B1[1]*C1[0];
+ tmp = sqrt(tmpVec * tmpVec);
+ for (int i=0; i<dim; ++i) { N1[i] = tmpVec[i] / tmp; }
+
+ // Calculate D1.
+ D1[0] = lambda*unit_B1[0] + nue*E_21[0] + dD*N1[0];
+ D1[1] = lambda*unit_B1[1] + nue*E_21[1] + dD*N1[1];
+ D1[2] = lambda*unit_B1[2] + nue*E_21[2] + dD*N1[2];
+
+// // Check projection of triangle A1B1C1 (-> triangle A2B2C2) to x-y-plane.
+// // Rotation is said to be not possible if
+// // -> angles are too small: | cos(A2B2,A2C2) | < 1-1.e-3
+// // -> area of A2B2C2 is too small:
+// // | det(A2B2 A2C2) | > max(1.e-10, 0.01*norm_B*norm_B)
+// double abs_det_B2C2 = fabs(B1[0]*C1[1]-B1[1]*C1[0]);
+// double limit = (0.01*norm_B*norm_B < 1.e-10) ?
+// 1.e-10 : 0.01*norm_B*norm_B;
+// if (abs_det_B2C2 < limit)
+// return false;
+
+// double norm_AB2 = sqrt(B1[0]*B1[0]+B1[1]*B1[1]);
+// double norm_AC2 = sqrt(C1[0]*C1[0]+C1[1]*C1[1]);
+// double cos_B2AC2 = (B1[0]*C1[0]+B1[1]*C1[1]) / (norm_AB2*norm_AC2);
+// if (fabs(cos_B2AC2) > 0.999)
+// return false;
+
+ return true;
+}
+
+const int
+ElementUpdate_3d::calcPosition(const WorldVector<double> &B2,
+ const WorldVector<double> &C2,
+ const WorldVector<double> &D2)
+{
+ FUNCNAME("ElementUpdate_3d::calcPosition");
+
+ double l_A, l_B2, l_C2;
+ double tmp_B2, tmp_C2;
+ double norm_AB2 = sqrt(B2 * B2);
+ double norm_AC2 = sqrt(C2 * C2);
+ double sP_AB2AC2 = B2 * C2;
+ double cos_B2AC2 = sP_AB2AC2 / (norm_AB2*norm_AC2);
+ int posFlag = -1;
+
+ // ===== Calculate barycentric coordinates of D2 with respect to
+ // triangle AB2C2. =====
+ if (fabs(cos_B2AC2) < 1.e-15) {
+ double sP_AB2AD2 = B2 * D2;
+ double sP_AC2AD2 = C2 * D2;
+
+ tmp_B2 = (sP_AB2AD2 < 0) ? -D2[0] : D2[0];
+ tmp_C2 = (sP_AC2AD2 < 0) ? -D2[1] : D2[1];
+
+ if (C2[1] < 0)
+ tmp_C2 *= -1;
+ }
+ else {
+ double sin_B2AC2 = sqrt(1-cos_B2AC2*cos_B2AC2);
+
+ if (C2[1] < 0)
+ sin_B2AC2 *= -1;
+
+ tmp_C2 = D2[1]/sin_B2AC2;
+ tmp_B2 = D2[0] - D2[1]*cos_B2AC2/sin_B2AC2;
+ }
+
+ l_B2 = tmp_B2 / norm_AB2;
+ l_C2 = tmp_C2 / norm_AC2;
+
+// // for test purposes: barycentric coordinates with Cramer's rule
+// double det_B2C2 = B2[0]*C2[1]-B2[1]*C2[0];
+
+// TEST_EXIT(fabs(det_B2C2) > 1.e-15)("illegal projection triangle !\n");
+
+// double det_D2C2 = D2[0]*C2[1]-D2[1]*C2[0];
+// double det_B2D2 = B2[0]*D2[1]-B2[1]*D2[0];
+
+// double l_B2_Cramer = 1.0*det_D2C2/det_B2C2;
+// double l_C2_Cramer = 1.0*det_B2D2/det_B2C2;
+
+// // TEST_EXIT(fabs(l_B2-l_B2_Cramer) < 1.e-8 &&
+// // fabs(l_C2-l_C2_Cramer) < 1.e-8)
+// TEST_EXIT(fabs(l_B2*l_B2_Cramer) >= 0 &&
+// fabs(l_C2*l_C2_Cramer) >= 0)
+// ("problems in barycentric coordinates of projection point !\n");
+// // end: for test purposes
+
+ /* // Corrections. */
+ /* if (fabs(l_B2) < 1.e-15) l_B2 = 0.0; */
+ /* if (fabs(1-l_B2) < 1.e-15) l_B2 = 1.0; */
+ /* if (fabs(l_C2) < 1.e-15) l_C2 = 0.0; */
+ /* if (fabs(1-l_C2) < 1.e-15) l_C2 = 1.0; */
+
+ l_A = 1-l_B2-l_C2;
+
+ //save barycentric coordinates for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setBarycentricCoords_3D(l_A,l_B2,l_C2,0);
+ }
+
+ // ===== Determine position of D2. =====
+ if (l_A < 0) {
+ if (l_B2 < 0) {
+ if (l_C2 < 0) {
+ ERROR_EXIT("calculation error: sum of barycentric coordinates is equal to 1 !\n");
+ }
+ else {
+ posFlag = VERT_C;
+ }
+ }
+ else if (l_C2 < 0) {
+ posFlag = VERT_B;
+ }
+ else { // EDGE_BC
+ posFlag = EDGE_BC;
+ /* if (l_B2 > 0) { */
+ /* posFlag = (l_C2 > 0) ? EDGE_BC : VERT_B; */
+ /* } */
+ /* else { */
+ /* if (l_C2 > 0) { */
+ /* posFlag = VERT_C; */
+ /* } */
+ /* else { */
+ /* ERROR_EXIT("no barycentric coordinates !\n"); */
+ /* } */
+ /* } */
+ }
+ }
+ else if (l_B2 < 0) {
+ if (l_C2 < 0) {
+ posFlag = VERT_A;
+ }
+ else { // EDGE_AC
+ posFlag = EDGE_AC;
+ /* if (l_A > 0) { */
+ /* posFlag = (l_C2 > 0) ? EDGE_AC : VERT_A; */
+ /* } */
+ /* else { */
+ /* if (l_C2 > 0) { */
+ /* posFlag = VERT_C; */
+ /* } */
+ /* else { */
+ /* ERROR_EXIT("no barycentric coordinates !\n"); */
+ /* } */
+ /* } */
+ }
+ }
+ else if (l_C2 < 0) { // EDGE_AB
+ posFlag = EDGE_AB;
+ /* if (l_A > 0) { */
+ /* posFlag = (l_B2 > 0) ? EDGE_AB : VERT_A; */
+ /* } */
+ /* else { */
+ /* if (l_B2 > 0) { */
+ /* posFlag = VERT_B; */
+ /* } */
+ /* else { */
+ /* ERROR_EXIT("no barycentric coordinates !\n"); */
+ /* } */
+ /* } */
+ }
+ else {
+ posFlag = IN_ABC;
+ }
+
+ return posFlag;
+}
+
+double
+ElementUpdate_3d::calcFaceUpdate(WorldVector<double> *A2d,
+ WorldVector<double> *B2d,
+ WorldVector<double> *C2d,
+ double &uhValA2d,
+ double &uhValB2d)
+{
+ int dim = A2d->getSize();
+ FixVec<WorldVector<double> *, VERTEX> vert2d(dim, NO_INIT);
+ FixVec<double, VERTEX> uhVal2d(dim, NO_INIT);
+
+ vert2d[0] = A2d;
+ vert2d[1] = B2d;
+ vert2d[2] = C2d;
+
+ uhVal2d[0] = uhValA2d;
+ uhVal2d[1] = uhValB2d;
+
+ return elUpdate2d->calcElementUpdate(vert2d, uhVal2d);
+}
diff --git a/AMDiS/Reinit/src/ElementUpdate_3d.h b/AMDiS/Reinit/src/ElementUpdate_3d.h
new file mode 100644
index 00000000..1801dd73
--- /dev/null
+++ b/AMDiS/Reinit/src/ElementUpdate_3d.h
@@ -0,0 +1,117 @@
+#ifndef ELEMENTUPDATE_3D_H
+#define ELEMENTUPDATE_3D_H
+
+#include "FixVec.h"
+// #include "MemoryManager.h"
+
+#include "ElementUpdate.h"
+#include "ElementUpdate_2d.h"
+#include "VelocityExt.h"
+
+using namespace AMDiS;
+
+class ElementUpdate_3d : public ElementUpdate
+{
+ public:
+ MEMORY_MANAGED(ElementUpdate_3d);
+
+ /**
+ * Constructor.
+ */
+ ElementUpdate_3d(VelocityExt *velExt_ = NULL)
+ : ElementUpdate(velExt_)
+ {
+ elUpdate2d = NEW ElementUpdate_2d(velExt_);
+ };
+
+ /**
+ * Destructor.
+ */
+ ~ElementUpdate_3d()
+ {
+ DELETE elUpdate2d;
+ };
+
+ /**
+ * Realization of ElementUpdate::calcElementUpdate.
+ * Calculates the Bornemann update for an element. The update is
+ * calculated for the last vertex of \ref vert, denoted by xh.
+ */
+ double calcElementUpdate(const FixVec<WorldVector<double> *, VERTEX> &vert_,
+ FixVec<double, VERTEX> &uhVal);
+
+ protected:
+ /**
+ * Sorts vertices of the face opposite to xh and translates the element
+ * to the point of origin.
+ *
+ * Denote A the vertex \ref vert[0], B the vertex \ref vert[1], ...
+ * after the sort, and uhVal(A), uhVal(B), ... the corresponding function
+ * values.
+ * Then, after the sort, uhVal(A) <= uhVal(B) <= uhVal(C).
+ * Finally, translation of element moves A to zero, uhVal(A) is set to
+ * zero, uhVal(B), uhVal(C) correspondingly.
+ */
+ void sortAndTranslateElement(FixVec<WorldVector<double> *, VERTEX> &vert,
+ FixVec<double, VERTEX> &uhVal,
+ double &uhVal0_orig);
+
+ /**
+ * Calculates rotation of element.
+ *
+ * Denote A1, B1, ... the vertices of element after rotation.
+ * Then distance of B1 to x-y-plane is uhVal(B) and
+ * distance of C1 to x-y-plane is uhVal(C).
+ *
+ * Return value: true - rotation exists
+ * false - rotation does not exist
+ *
+ * Assumptions: --> A == zero vector (then AB == B, ...)
+ * --> uhVal(A) == 0
+ */
+ bool rotateElement(const FixVec<WorldVector<double> *, VERTEX> &vert,
+ const FixVec<double, VERTEX> &uhVal,
+ WorldVector<double> &B1,
+ WorldVector<double> &C1,
+ WorldVector<double> &D1);
+
+ /**
+ * Determines position of D2 with respect to the triangle AB2C2.
+ * A is the point of origin.
+ *
+ * Note: Limited to triangles with angles not greater than Pi/2
+ * (acute triangles).
+ */
+ const int calcPosition(const WorldVector<double> &B2,
+ const WorldVector<double> &C2,
+ const WorldVector<double> &D2);
+
+ /**
+ * Calculate 2-dimensional Bornemann update for face of element. The
+ * update is calculated for vertex with coordinates \ref *C2d.
+ */
+ double calcFaceUpdate(WorldVector<double> *A2d,
+ WorldVector<double> *B2d,
+ WorldVector<double> *C2d,
+ double &uhValA2d,
+ double &uhValB2d);
+
+ protected:
+ /**
+ * Used for Bornemann update on faces of element.
+ */
+ ElementUpdate_2d *elUpdate2d;
+
+ /**
+ * Flags to decribe position of D2 with respect to AB2C2.
+ */
+ static const int IN_ABC = 0;
+ static const int VERT_A = 1;
+ static const int VERT_B = 2;
+ static const int VERT_C = 3;
+ static const int EDGE_AB = 4;
+ static const int EDGE_AC = 5;
+ static const int EDGE_BC = 6;
+};
+
+#endif // ELEMENTUPDATE_3D_H
diff --git a/AMDiS/Reinit/src/HL_SignedDist.cc b/AMDiS/Reinit/src/HL_SignedDist.cc
new file mode 100644
index 00000000..a81ae35a
--- /dev/null
+++ b/AMDiS/Reinit/src/HL_SignedDist.cc
@@ -0,0 +1,375 @@
+#include "HL_SignedDist.h"
+
+const Flag HL_SignedDist::VEL_EXT = 0X01L;
+const Flag HL_SignedDist::VEL_EXT_FROM_VEL_FIELD = 0X02L;
+
+void
+HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ ElementFunction<double> *elFct)
+{
+ adaptInfo = adaptInfo_;
+
+ if (elFct == NULL) {
+ TEST_EXIT(lS_DOF_)("illegal level set function !\n");
+ TEST_EXIT(lS_DOF_->getFESpace() == sD_DOF_->getFESpace())("DOF vectors do not have the same FE space!\n");
+ }
+ TEST_EXIT(sD_DOF_)("illegal DOF vector for signed distance function !\n");
+ TEST_EXIT(sD_DOF_->getFESpace()->getBasisFcts()->getDegree() == 1)("does it work for non-linear finite elements ?\n");
+
+ lS_DOF = lS_DOF_;
+ sD_DOF = sD_DOF_;
+ feSpace = sD_DOF->getFESpace();
+
+ // ===== Initialization. =====
+ initialize(elFct);
+
+ // ===== Print level set function. =====
+ if (lS_DOF)
+ printLevelSetFct();
+
+ // ===== Create DOF vector to mark boundary vertices. =====
+ if (bound_DOF)
+ DELETE bound_DOF;
+ bound_DOF = NEW DOFVector<double>(feSpace, "bound_DOF");
+ bound_DOF->set(0.0);
+
+ // ===== Boundary vertex and boundary value initialization. =====
+ initializeBoundary();
+
+ // ===== Print boundary initialization. =====
+ printBoundInitFct();
+
+ // ===== Calculate distance function. =====
+ HL_updateIteration();
+
+ // ===== Transformation to signed distance function. =====
+ if (elFct == NULL)
+ setSign();
+
+ // ===== Print calculated signed distance function. =====
+ printSignedDistFct();
+}
+
+void
+HL_SignedDist::calcSignedDistFct(AdaptInfo *adaptInfo_,
+ DOFVector<double> *lS_DOF_)
+{
+ TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
+
+ sD_DOF = NEW DOFVector<double>(lS_DOF_->getFESpace(), "sD_DOF");
+
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, NULL);
+
+ // lS_DOF_->copy(*sD_DOF);
+ *lS_DOF_ = *sD_DOF;
+
+ DELETE sD_DOF;
+}
+
+void
+HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
+ DOFVector<double> *origVel_DOF_,
+ DOFVector<double> *vel_DOF_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct)
+{
+ TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
+ TEST_EXIT(velExt)("velExt not defined !\n");
+
+ velExt->setVelocity(origVel_DOF_, vel_DOF_);
+
+ velExt->printOrigVelDOF(adaptInfo_);
+
+ if (calcSDFct || sD_DOF_ != NULL) {
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
+ }
+ else {
+
+ TEST_EXIT(vel_DOF_)("illegal level set function lS_DOF_ !\n");
+
+ sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
+
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
+
+ DELETE sD_DOF;
+ }
+
+ velExt->printVelDOF(adaptInfo_);
+}
+
+void
+HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
+ DOFVector<double> *vel_DOF_,
+ DOFVector<double> *lS_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct)
+{
+ TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
+ TEST_EXIT(velExt)("velExt not defined !\n");
+ TEST_EXIT(vel_DOF_)("velocity vector vel_DOF_ not defined !\n");
+
+ if (calcSDFct)
+ TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
+
+ DOFVector<double> *newVel_DOF_ =
+ NEW DOFVector<double>(vel_DOF_->getFESpace(), "vel_DOF_");
+
+ velExt->setVelocity(vel_DOF_, newVel_DOF_);
+
+ velExt->printOrigVelDOF(adaptInfo_);
+
+ if (calcSDFct && elFct == NULL)
+ calcSignedDistFct(adaptInfo_, lS_DOF_);
+ else {
+
+ sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
+
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
+
+ if (calcSDFct)
+ lS_DOF_->copy(*sD_DOF);
+
+ DELETE sD_DOF;
+ }
+
+ vel_DOF_->copy(*newVel_DOF_);
+
+ velExt->printVelDOF(adaptInfo_);
+
+ DELETE newVel_DOF_;
+}
+
+void
+HL_SignedDist::calcVelocityExt(
+ AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> origVel_DOF_,
+ std::vector<DOFVector<double> *> vel_DOF_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct)
+{
+ TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
+ TEST_EXIT(velExt)("velExt not defined !\n");
+
+ velExt->setVelocity(origVel_DOF_, vel_DOF_);
+
+ velExt->printOrigVelDOF(adaptInfo_);
+
+ if (calcSDFct || sD_DOF_ != NULL) {
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
+ }
+ else {
+
+ TEST_EXIT(vel_DOF_[0])("illegal level set function lS_DOF_ !\n");
+
+ sD_DOF = NEW DOFVector<double>(vel_DOF_[0]->getFESpace(), "sD_DOF");
+
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
+
+ DELETE sD_DOF;
+ }
+
+ velExt->printVelDOF(adaptInfo_);
+}
+
+void
+HL_SignedDist::calcVelocityExt(AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> vel_DOF_,
+ DOFVector<double> *lS_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct)
+{
+ TEST_EXIT(velExtType.isSet(VEL_EXT))("illegal velocity extension type !\n");
+ TEST_EXIT(velExt)("velExt not defined !\n");
+
+ int nVelDOFs = vel_DOF_.size();
+
+ if (calcSDFct)
+ TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
+
+ std::vector<DOFVector<double> *> newVel_DOF_(nVelDOFs);
+ for (int nV=0; nV<nVelDOFs; ++nV) {
+
+ TEST_EXIT(vel_DOF_[nV])("velocity vector vel_DOF_ not defined !\n");
+
+ newVel_DOF_[nV] =
+ NEW DOFVector<double>((vel_DOF_[nV])->getFESpace(), "vel_DOF_");
+ }
+
+ velExt->setVelocity(vel_DOF_, newVel_DOF_);
+
+ velExt->printOrigVelDOF(adaptInfo_);
+
+ if (calcSDFct && elFct == NULL)
+ calcSignedDistFct(adaptInfo_, lS_DOF_);
+ else {
+
+ sD_DOF = NEW DOFVector<double>(vel_DOF_[0]->getFESpace(), "sD_DOF");
+
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
+
+ if (calcSDFct)
+ lS_DOF_->copy(*sD_DOF);
+
+ DELETE sD_DOF;
+ }
+
+ for (int nV=0; nV<nVelDOFs; ++nV)
+ (vel_DOF_[nV])->copy(*(newVel_DOF_[nV]));
+
+ velExt->printVelDOF(adaptInfo_);
+
+ for (int nV=0; nV<nVelDOFs; ++nV)
+ DELETE newVel_DOF_[nV];
+}
+
+void
+HL_SignedDist::calcVelocityExtFromVelocityField(
+ AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> &velField_,
+ DOFVector<double> *vel_DOF_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct)
+{
+ TEST_EXIT(velExtType.isSet(VEL_EXT_FROM_VEL_FIELD))
+ ("illegal velocity extension type !\n");
+ TEST_EXIT(velExt)("velExt not defined !\n");
+ TEST_EXIT(elFct == NULL)("not implemented yet for elFct != NULL !\n");
+
+ ((VelocityExtFromVelocityField *)(velExt))->setVelocityField(velField_, lS_DOF_, vel_DOF_);
+
+ if (calcSDFct) {
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF_, elFct);
+ }
+ else {
+
+ TEST_EXIT(vel_DOF_)("illegal velocity vector vel_DOF_ !\n");
+
+ sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
+
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
+
+ DELETE sD_DOF;
+ }
+
+ velExt->printVelDOF(adaptInfo_);
+}
+
+void
+HL_SignedDist::calcVelocityExtFromVelocityField(
+ AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> &velField_,
+ DOFVector<double> *vel_DOF_,
+ DOFVector<double> *lS_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct)
+{
+ TEST_EXIT(velExtType.isSet(VEL_EXT_FROM_VEL_FIELD))
+ ("illegal velocity extension type !\n");
+ TEST_EXIT(velExt)("velExt not defined !\n");
+ TEST_EXIT(elFct == NULL)("not implemented yet for elFct != NULL !\n");
+
+ if (calcSDFct)
+ TEST_EXIT(lS_DOF_)("illegal level set function lS_DOF_ !\n");
+
+ ((VelocityExtFromVelocityField *)(velExt))->setVelocityField(velField_, lS_DOF_, vel_DOF_);
+
+ if (calcSDFct && elFct == NULL) {
+ calcSignedDistFct(adaptInfo_, lS_DOF_);
+ }
+ else {
+
+ TEST_EXIT(vel_DOF_)("illegal velocity vector vel_DOF_ !\n");
+
+ sD_DOF = NEW DOFVector<double>(vel_DOF_->getFESpace(), "sD_DOF");
+
+ calcSignedDistFct(adaptInfo_, lS_DOF_, sD_DOF, elFct);
+
+ if (calcSDFct)
+ lS_DOF_->copy(*sD_DOF);
+
+ DELETE sD_DOF;
+ }
+
+ velExt->printVelDOF(adaptInfo_);
+}
+
+void
+HL_SignedDist::initialize(ElementFunction<double> *elFct)
+{
+ if (elUpdate)
+ DELETE elUpdate;
+ if (bndElDist)
+ DELETE bndElDist;
+
+ if (elLS)
+ DELETE elLS;
+ if (phi)
+ DELETE phi;
+
+ // ===== Create ElementLevelSet. =====
+ if (elFct == NULL) {
+ phi = NEW ElementFunctionDOFVec<double>(lS_DOF);
+ elLS = NEW ElementLevelSet("ElementLevelSet", phi, feSpace->getMesh());
+ }
+ else {
+ elLS = NEW ElementLevelSet("ElementLevelSet", elFct, feSpace->getMesh());
+ }
+
+
+ // ===== Define boundary initialization strategy. =====
+ int boundInitFlag;
+ GET_PARAMETER(0, name + "->boundary initialization", "%d",
+ &boundInitFlag);
+ switch(boundInitFlag) {
+ case 0: bndElDist = NEW BoundaryElementLevelSetDist(elLS, dim);
+ break;
+ case 1: bndElDist = NEW BoundaryElementTopDist(elLS, dim, velExt);
+ break;
+ case 2: bndElDist = NEW BoundaryElementEdgeDist(elLS, dim);
+ break;
+ case 3: bndElDist = NEW BoundaryElementNormalDist(elLS, dim);
+ break;
+ default: ERROR_EXIT("illegal boundary initialization !\n");
+ break;
+ }
+
+ // ===== Create elUpdate. =====
+ TEST_EXIT(velExt == NULL || boundInitFlag == 1)
+ ("velocity extension only works with topological boundary element initialization method !\n");
+
+ switch(dim) {
+ case 2: elUpdate = NEW ElementUpdate_2d(velExt);
+ break;
+ case 3: elUpdate = NEW ElementUpdate_3d(velExt);
+ break;
+ default: ERROR_EXIT("illegal dimension !\n");
+ break;
+ }
+}
+
+void
+HL_SignedDist::setSign()
+{
+ DOFVector<double>::Iterator it_sD(const_cast<DOFVector<double> *>(sD_DOF),
+ USED_DOFS);
+ DOFVector<double>::Iterator it_lS(const_cast<DOFVector<double> *>(lS_DOF),
+ USED_DOFS);
+
+ for(it_sD.reset(), it_lS.reset();
+ !it_sD.end();
+ ++it_sD, ++it_lS) {
+
+ // Vertex lies inside the zero level set ?
+ if ((*it_lS) < 0) {
+ (*it_sD) *= -1;
+ }
+ }
+}
diff --git a/AMDiS/Reinit/src/HL_SignedDist.h b/AMDiS/Reinit/src/HL_SignedDist.h
new file mode 100644
index 00000000..0853c6f4
--- /dev/null
+++ b/AMDiS/Reinit/src/HL_SignedDist.h
@@ -0,0 +1,407 @@
+#ifndef HL_SIGNEDDIST
+#define HL_SIGNEDDIST
+
+#include "Global.h"
+#include "AdaptInfo.h"
+#include "DOFVector.h"
+#include "ElementFunction.h"
+#include "FileWriter.h"
+#include "FixVec.h"
+#include "Flag.h"
+// #include "MemoryManager.h"
+#include "Parameters.h"
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+#include "BoundaryElementLevelSetDist.h"
+#include "BoundaryElementTopDist.h"
+#include "BoundaryElementEdgeDist.h"
+#include "BoundaryElementNormalDist.h"
+#include "ElementUpdate.h"
+#include "ElementUpdate_2d.h"
+#include "ElementUpdate_3d.h"
+#include "VelocityExt.h"
+#include "VelocityExtFromVelocityField.h"
+
+using namespace AMDiS;
+
+//////////////////////////////////////////////////////////////////////////////
+//
+// class HL_SignedDist:
+// --------------------
+// Holds functionality for
+// - calculation of signed distance function for a level set function
+// (calcSignedDistFct())
+// - extension of velocity from an interface on complete mesh in direction
+// normal to the interface, where the interface is given implicitly
+// by a level set function
+// (calcVelocityExt())
+//
+//////////////////////////////////////////////////////////////////////////////
+class HL_SignedDist
+{
+ public:
+ MEMORY_MANAGED(HL_SignedDist);
+
+ /**
+ * Constructor.
+ */
+ HL_SignedDist(const char *name_,
+ int dim_,
+ bool doVelocityExt = false,
+ Flag velExtType_ = VEL_EXT)
+ : name(name_),
+ adaptInfo(NULL),
+ dim(dim_),
+ lS_DOF(NULL),
+ sD_DOF(NULL),
+ bound_DOF(NULL),
+ feSpace(NULL),
+ elUpdate(NULL),
+ bndElDist(NULL),
+ elLS(NULL),
+ phi(NULL),
+ velExt(NULL),
+ velExtType(velExtType_)
+ {
+ FUNCNAME("HL_SignedDist::HL_SignedDist");
+
+ TEST_EXIT(dim == 2 || dim == 3)("only works for dimension 2 and 3 !\n");
+
+ // ===== Read parameters from init file. =====
+ GET_PARAMETER(0, name + "->infinity value", "%f", &inftyValue);
+
+ TEST_EXIT(inftyValue > 1000)("illegal infinity value !\n");
+
+ // ===== Create functionality for velocity extension. =====
+ if (doVelocityExt) {
+ if (velExtType.isSet(VEL_EXT))
+ velExt = NEW VelocityExt(dim);
+ else
+ velExt = NEW VelocityExtFromVelocityField(dim);
+ }
+ };
+
+ /**
+ * Virtual destructor.
+ */
+ virtual ~HL_SignedDist()
+ {
+ if (elUpdate)
+ DELETE elUpdate;
+ if (bndElDist)
+ DELETE bndElDist;
+
+ if (elLS)
+ DELETE elLS;
+ if (phi)
+ DELETE phi;
+
+ if (bound_DOF)
+ DELETE bound_DOF;
+
+ DELETE velExt;
+ };
+
+ /**
+ * Calculates the signed distance function for the interface given
+ * implicitly by the zero level set of lS_DOF_. The result
+ * is stored in sD_DOF_.
+ *
+ * Alternative to level set function as DOF vector:
+ * If elFct != NULL, this ElementFunction is used as level set function.
+ * In this case: only calculation of distance function (positive sign) !
+ */
+ void calcSignedDistFct(AdaptInfo *adaptInfo_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Calculates the signed distance function for the interface given
+ * implicitly by the zero level set of lS_DOF_. The result
+ * is stored in lS_DOF_.
+ */
+ void calcSignedDistFct(AdaptInfo *adaptInfo_,
+ DOFVector<double> *lS_DOF_);
+
+ /**
+ * Calculates the extension of a velocity origVel_DOF_ from the interface
+ * on the complete mesh in direction normal to the interface. The interface
+ * is given implicitly as the zero level set of lS_DOF_. The result
+ * is stored in vel_DOF_. If calcSDFct is true, the signed distance function
+ * which is calculated during the extension of the velocity anyway is
+ * stored in sD_DOF_.
+ *
+ * Alternative to level set function as DOF vector:
+ * If elFct != NULL, this ElementFunction is used as level set function.
+ * In this case: only calculation of distance function (positive sign) !
+ */
+ void calcVelocityExt(AdaptInfo *adaptInfo_,
+ DOFVector<double> *origVel_DOF_,
+ DOFVector<double> *vel_DOF_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Calculates the extension of a velocity vel_DOF_ from the interface
+ * on the complete mesh in direction normal to the interface. The interface
+ * is given implicitly as the zero level set of lS_DOF_. The result
+ * is stored in origVel_DOF_. If calcSDFct is true, the signed distance
+ * function, which is calculated during the extension of the velocity
+ * anyway, is stored in lS_DOF_.
+ *
+ * Alternative to level set function as DOF vector:
+ * If elFct != NULL, this ElementFunction is used as level set function.
+ * In this case: only calculation of distance function (positive sign) !
+ */
+ void calcVelocityExt(AdaptInfo *adaptInfo_,
+ DOFVector<double> *vel_DOF_,
+ DOFVector<double> *lS_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Calculates the extension of the velocity vectors in origVel_DOF_
+ * from the interface on the complete mesh in direction normal to the
+ * interface. The interface is given implicitly as the zero level set of
+ * lS_DOF_. The result is stored in vel_DOF_.
+ * If calcSDFct is true, the signed distance function,
+ * which is calculated during the extension of the velocity anyway, is
+ * stored in sD_DOF_.
+ *
+ * Alternative to level set function as DOF vector:
+ * If elFct != NULL, this ElementFunction is used as level set function.
+ * In this case: only calculation of distance function (positive sign) !
+ */
+ void calcVelocityExt(AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> origVel_DOF_,
+ std::vector<DOFVector<double> *> vel_DOF_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Calculates the extension of the velocity vectors in vel_DOF_
+ * from the interface on the complete mesh in direction normal to the
+ * interface. The interface is given implicitly as the zero level set
+ * of lS_DOF_. The result is stored in vel_DOF_.
+ * If calcSDFct is true, the signed distance function, which is calculated
+ * during the extension of the velocity anyway, is stored in lS_DOF_.
+ *
+ * Alternative to level set function as DOF vector:
+ * If elFct != NULL, this ElementFunction is used as level set function.
+ * In this case: only calculation of distance function (positive sign) !
+ */
+ void calcVelocityExt(AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> vel_DOF_,
+ DOFVector<double> *lS_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Calculates the normal velocity for the velocity field velField_ and its
+ * extension from the interface on the complete mesh in direction normal
+ * to the interface. The interface is given implicitly as the
+ * zero level set of lS_DOF_. The result is stored in vel_DOF_.
+ * If calcSDFct is true, the signed distance function
+ * which is calculated during the extension of the velocity anyway is
+ * stored in sD_DOF_.
+ *
+ * Alternative to level set function as DOF vector:
+ * If elFct != NULL, this ElementFunction is used as level set function.
+ * In this case: only calculation of distance function (positive sign) !
+ */
+ void calcVelocityExtFromVelocityField(
+ AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> &velField_,
+ DOFVector<double> *vel_DOF_,
+ const DOFVector<double> *lS_DOF_,
+ DOFVector<double> *sD_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Calculates the normal velocity for the velocity field velField_ and its
+ * extension from the interface on the complete mesh in direction normal
+ * to the interface. The interface is given implicitly as the
+ * zero level set of lS_DOF_. The result is stored in vel_DOF_.
+ * If calcSDFct is true, the signed distance function
+ * which is calculated during the extension of the velocity anyway is
+ * stored in lS_DOF_.
+ *
+ * Alternative to level set function as DOF vector:
+ * If elFct != NULL, this ElementFunction is used as level set function.
+ * In this case: only calculation of distance function (positive sign) !
+ */
+ void calcVelocityExtFromVelocityField(
+ AdaptInfo *adaptInfo_,
+ std::vector<DOFVector<double> *> &velField_,
+ DOFVector<double> *vel_DOF_,
+ DOFVector<double> *lS_DOF_,
+ bool calcSDFct,
+ ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Print initial function: level set function defining the interface.
+ */
+ void printLevelSetFct()
+ {
+ FileWriter *fileWriter = NEW FileWriter(
+ "SignedDist->level set fct output",
+ feSpace->getMesh(),
+ const_cast<DOFVector<double> *>(lS_DOF));
+ fileWriter->writeFiles(adaptInfo, false);
+
+ DELETE fileWriter;
+ };
+
+ /**
+ * Print signed distance function.
+ */
+ void printSignedDistFct()
+ {
+ FileWriter *fileWriter = NEW FileWriter(
+ "SignedDist->result output",
+ feSpace->getMesh(),
+ sD_DOF);
+ fileWriter->writeFiles(adaptInfo, false);
+
+ DELETE fileWriter;
+ };
+
+ protected:
+ /**
+ * Initialization.
+ */
+ virtual void initialize(ElementFunction<double> *elFct = NULL);
+
+ /**
+ * Initializes the boundary: calculation of the distance of boundary
+ * vertices to the interface.
+ * Interface is given by lS_DOF and result is stored in
+ * sD_DOF.
+ */
+ virtual void initializeBoundary() = 0;
+
+ /**
+ * Calculates the distance function and stores result in sD_DOF.
+ * Requirement: The boundary values are already set in sD_DOF.
+ */
+ virtual void HL_updateIteration() = 0;
+
+ /**
+ * Transforms the distance function into a signed distance function.
+ * The sign is given by the level set function lS_DOF. The
+ * signed distance function is stored in sD_DOF.
+ */
+ void setSign();
+
+ /**
+ * Print boundary initialization (initial function for Hopf-Lax iteration).
+ */
+ void printBoundInitFct()
+ {
+ FileWriter *fileWriter = NEW FileWriter(
+ "SignedDist->boundary initialization output",
+ feSpace->getMesh(),
+ sD_DOF);
+ fileWriter->writeFiles(adaptInfo, false);
+
+ DELETE fileWriter;
+ };
+
+ public:
+ /**
+ * Flags to distinguish velocity extension types.
+ */
+ static const Flag VEL_EXT;
+ static const Flag VEL_EXT_FROM_VEL_FIELD;
+
+ protected:
+ /**
+ * Name of this instantiation of HL_SignedDist.
+ */
+ std::string name;
+
+ /**
+ * AdaptInfo.
+ */
+ AdaptInfo *adaptInfo;
+
+ /**
+ * Dimension.
+ */
+ int dim;
+
+ /**
+ * Level set function giving implicitely (as zero level set) the
+ * interface for which the signed distance function is calculated.
+ */
+ const DOFVector<double> *lS_DOF;
+
+ /**
+ * DOF vector for the calculated signed distance function.
+ * Also used during calculation.
+ */
+ DOFVector<double> *sD_DOF;
+
+ /**
+ * Marker for boundary vertices:
+ * 0 - vertex is no boundary vertex
+ * 1 - vertex is boundary vertex
+ */
+ DOFVector<double> *bound_DOF;
+
+ /**
+ * Finite element space.
+ */
+ const FiniteElemSpace *feSpace;
+
+ /**
+ * Initialization value "inifinity" for non-boundary vertices.
+ */
+ double inftyValue;
+
+ /**
+ * Pointer to ElementUpdate. Used for Hopf-Lax element update.
+ */
+ ElementUpdate *elUpdate;
+
+ /**
+ * Used for boundary vertex initialization: calculation of the distance
+ * to the interface for all vertices of a boundary element.
+ */
+ BoundaryElementDist *bndElDist;
+
+ /**
+ * Holds level set function and functionalities for intersection point
+ * calculation.
+ */
+ ElementLevelSet *elLS;
+
+ /*
+ * Level set function which implicitely gives the interface as its
+ * zero level set.
+ * This representation is needed for the use of class ElementLevelSet.
+ */
+ ElementFunction<double> *phi;
+
+ /**
+ * Object needed to extrapolate velocity from the interface.
+ */
+ VelocityExt *velExt;
+
+ /**
+ * Type of velocity extension method. Possible types:
+ * VEL_EXT - object of class VelocityExt
+ * VEL_EXT_FROM_VEL_FIELD - object of class VelocityExtFromVelocityField
+ */
+ Flag velExtType;
+};
+
+#endif // HL_SIGNEDDIST
diff --git a/AMDiS/Reinit/src/HL_SignedDistBornemann.h b/AMDiS/Reinit/src/HL_SignedDistBornemann.h
new file mode 100644
index 00000000..648e97a7
--- /dev/null
+++ b/AMDiS/Reinit/src/HL_SignedDistBornemann.h
@@ -0,0 +1,916 @@
+#ifndef HL_SIGNEDDISTBORNEMANN
+#define HL_SIGNEDDISTBORNEMANN
+
+#include "ElInfo.h"
+#include "FixVec.h"
+// #include "MemoryManager.h"
+#include "Traverse.h"
+#include <queue>
+#include <time.h>
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+/* #include "BoundaryElementTopDist.h" */
+/* #include "BoundaryElementEdgeDist.h" */
+/* #include "BoundaryElementNormalDist.h" */
+#include "ElementUpdate.h"
+#include "ElementUpdate_2d.h"
+#include "ElementUpdate_3d.h"
+#include "HL_SignedDist.h"
+
+#include "SMIAdapter.h"
+#include "smi.h"
+
+using namespace AMDiS;
+
+typedef struct
+{
+ int ElNum;
+ int VertNum;
+} Vert_Struct;
+
+class HL_SignedDistBornemann : public HL_SignedDist
+{
+ public:
+ MEMORY_MANAGED(HL_SignedDistBornemann);
+
+ /**
+ * Constructor.
+ */
+ HL_SignedDistBornemann(const char *name_,int dim_)
+ : HL_SignedDist(name_, dim_),
+ smiAdapter(NULL)
+ {
+ FUNCNAME("HL_SignedDistBornemann::HL_SignedDistBornemann");
+
+ // ===== Read parameters from init file. =====
+ GET_PARAMETER(0,name + "->tolerance", "%f", &tol);
+ GET_PARAMETER(0,name + "->count_how_often_saved_in_list", "%d", &count_in_list);
+ GET_PARAMETER(0,name + "->save_in_list->the ..th", "%d", &print_in_list);
+
+ TEST_EXIT(tol > 0)("illegal tolerance !\n");
+ };
+
+ protected:
+ /**
+ * Initializes the boundary: calculation of the distance of boundary
+ * vertices to the interface.
+ * Interface is given by \ref lS_DOF and result is stored in
+ * \ref sD_DOF.
+ */
+ void initializeBoundary()
+ {
+ FUNCNAME("HL_SignedDistBornemann::initializeBoundary()");
+
+ TraverseStack stack;
+ FixVec<double, VERTEX> distVec(dim, NO_INIT);
+ int elStatus;
+ const int *elVertStatusVec;
+ int numIntersecPoints;
+ int NumVertIntPoints;
+ int ElNum;
+ clock_t time1;
+ clock_t time2;
+ double TIME;
+
+ //=== transvering Mesh to SMI and add quantities
+ time1 = clock();
+ Mesh_to_SMI_and_quantity ();
+ time2 = clock();
+ TIME = TIME_USED(time1, time2);
+
+ cout <<"Zeit zum Transformieren nach SMI: "<<TIME<<" sec.\n";
+
+ // ===== All non-boundary vertices are initialized with "infinity". =====
+ sD_DOF->set(inftyValue);
+
+ // ===== Traverse mesh and initialize boundary elements. =====
+ const int nBasFcts = feSpace->getBasisFcts()->getNumber();
+ DegreeOfFreedom *locInd = GET_MEMORY(DegreeOfFreedom, nBasFcts);
+
+ ElInfo *elInfo = stack.traverseFirst(feSpace->getMesh(),
+ -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_COORDS);
+ while(elInfo) {
+
+ // Get local indices of vertices.
+ feSpace->getBasisFcts()->getLocalIndices(
+ const_cast<Element *>(elInfo->getElement()),
+ const_cast<DOFAdmin *>(feSpace->getAdmin()),
+ locInd);
+
+ // Get element status.
+ elStatus = ElementLevelSet::createElementLevelSet(elInfo);
+
+ //Get some information for creating the first list
+
+ elVertStatusVec = ElementLevelSet::getElVertStatusVec();
+ numIntersecPoints = ElementLevelSet::getNumElIntersecPoints();
+
+ //Beginn creating the first list
+
+ NumVertIntPoints = ElementLevelSet::getNumVertIntPoints();
+ ElNum = elInfo->getElement()->getIndex();
+
+ createLists( elStatus, elVertStatusVec, numIntersecPoints, ElNum, NumVertIntPoints, locInd);
+
+ //end creating two lists
+
+ // Is element cut by the interface ?
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+
+ // Reset element distance vector.
+ for (int i=0; i<=dim; ++i) {
+ distVec[i] = inftyValue;
+ }
+
+ // Mark all vertices as boundary vertices.
+ for (int i=0; i<=dim; ++i) {
+ (*bound_DOF)[locInd[i]] = 1.0;
+ }
+
+
+
+ // Calculate distance for all vertices.
+ if (bndElDist->calcDistOnBoundaryElement(elInfo, distVec) !=
+ ElementLevelSet::LEVEL_SET_BOUNDARY) {
+ ERROR_EXIT("error in distance calculation !\n");
+ }
+ else {
+
+ // If distance is smaller, correct to new distance.
+ for (int i=0; i<=dim; ++i) {
+ if ((*sD_DOF)[locInd[i]] > distVec[i]) {
+ (*sD_DOF)[locInd[i]] = distVec[i];
+ }
+ }
+ }
+ } // end of: elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY
+
+ else if (ElementLevelSet::getNumVertIntPoints() != 0){
+ // Interface cuts element only in vertices.
+ elVertStatusVec = ElementLevelSet::getElVertStatusVec();
+
+ for (int i=0; i<=dim; ++i) {
+ if (elVertStatusVec[i] == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+ (*bound_DOF)[locInd[i]] = 1.0;
+ (*sD_DOF)[locInd[i]] = 0.0;
+ }
+ }
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ } // end of: mesh traverse
+
+ FREE_MEMORY(locInd, DegreeOfFreedom, nBasFcts);
+
+ return;
+ };
+
+ /**
+ * Calculates the distance function and stores result in sD_DOF.
+ * Requirement: The boundary values are already set in \ref sD_DOF.
+ */
+ void HL_updateIteration()
+ {
+ int numNodes;
+ int *nodeIndices;
+ int max_q3 = 0;
+ int max_q4 = 0;
+ int sum_q3 = 0;
+ int sum_q4 = 0;
+ clock_t time1;
+ clock_t time2;
+ double TIME;
+
+ SMI_Get_all_nodes(1,1,const_cast<DegreeOfFreedom*>(&numNodes),&nodeIndices);
+ int *values_q3 = new int [numNodes];
+ int *values_q4 = new int [numNodes];
+
+ GET_PARAMETER(0,"SignedDist->count updates", "%d", &count_updates);
+
+ time1 = clock();
+ traverseListElement();
+ time2 = clock();
+ TIME = TIME_USED(time1, time2);
+
+ cout<<"Zeit zum durchlaufen der ersten Liste: "<<TIME<<" sec.\n\n";
+
+ time1 = clock();
+ traversingListELVert ( sD_DOF );
+ time2 = clock();
+ TIME = TIME_USED(time1,time2);
+
+ cout<<"Zeit zum Durchlaufen der zweiten Liste: "<<TIME<<" sec.\n\n";
+
+ std::string smiOutFile;
+ GET_PARAMETER(0, "SignedDist->count updates->output->filename", &smiOutFile);
+ cout << "count updates Ausgabe-Datei: " << smiOutFile.c_str() << "\n\n";
+
+ ofstream out (smiOutFile.c_str());
+
+ SMI_Get_quantity_values(1,1,3,SMI_TYPE_INT,1,numNodes,nodeIndices,values_q3);
+ SMI_Get_quantity_values(1,1,4,SMI_TYPE_INT,1,numNodes,nodeIndices,values_q4);
+
+ for (int i=0; i<numNodes; i++)
+ {
+ out<<nodeIndices[i]<<" "<<values_q3[i]<<" "<<values_q4[i]<<"\n";
+
+ if(max_q3 < values_q3[i])
+ {
+ max_q3 = values_q3[i];
+ }
+ if(max_q4 < values_q4[i])
+ {
+ max_q4 = values_q4[i];
+ }
+ sum_q3 = sum_q3 + values_q3[i];
+ sum_q4 = sum_q4 + values_q4[i];
+ }
+
+ out<<"\n\n maximale Anzahl an versuchten Updates auf einem Knoten: "<<max_q3<<"\n maximale Anzahl an durchgefuehrten Updates auf einem Knoten: "<<max_q4<<"\n";
+ out<<"\n Summe aller versuchten Updates: "<<sum_q3<<"\n Summe aller durchgefuehrten updates: "<<sum_q4<<"\n";
+
+ out<<"Anzahl an Knoten: "<<numNodes<<"\n";
+
+ out.close();
+
+
+
+
+ SMI_End_write_transaction(1,1);
+
+ return;
+ };
+
+ /**
+ *function for transvering Mesh to SMi an Adding the quantities
+ **/
+
+ void Mesh_to_SMI_and_quantity ()
+ {
+ if(!smiAdapter)
+ {
+ smiAdapter = NEW SMIAdapter(1,1,const_cast<FiniteElemSpace*>(feSpace),1,-1);
+ }
+
+ // cout << "\n\n\tSMI-Adapter angelegt !\n\n";
+
+ //====== transfer Mesh to SMI ======================
+
+ smiAdapter->addNeighbourInfo();
+
+ // cout << "\n\n\tNachbarschafts-Infos in SMI ergaenzt !\n\n";
+
+ smiAdapter-> transferMeshToSMI();
+
+ // cout << "\n\n\tGitter nach SMI geschrieben !\n\n";
+
+ SMI_Begin_write_transaction(1,1);
+
+ int nul=0;
+ //which pair of element and node is saved in list
+ SMI_Add_quantity(1, 1, 1, SMI_LOC_ELEM, SMI_TYPE_INT,dim+1, &nul);
+ //which node is a boundary-node
+ SMI_Add_quantity(1,1,2,SMI_LOC_NODE, SMI_TYPE_INT,1,&nul);
+ //saves the number of tried updates
+ SMI_Add_quantity(1,1,3,SMI_LOC_NODE,SMI_TYPE_INT,1,&nul);
+ //saves the number of updates
+ SMI_Add_quantity(1,1,4,SMI_LOC_NODE,SMI_TYPE_INT,1,&nul);
+ //how often a node is saved in the second list
+ SMI_Add_quantity(1,1,5,SMI_LOC_NODE,SMI_TYPE_INT,1,&nul);
+
+ // cout << "\n\n\tQuantities in SMI ergaenzt !\n\n";
+ }
+
+ //Begin fuction for creating two lists
+ //*
+
+ void createLists(const int elStatus,
+ const int *elVertStatusVec,
+ const int numIntersecPoints,
+ const int ElNum,
+ const int NumVertIntPoints,
+ const DegreeOfFreedom *locInd)
+ {
+ Vert_Struct vertStruct;
+ int *sv = new int [dim+1];
+ int value=0;
+ int value_q5 = 0;
+ int * nodeIndices;
+ /* int *nodes = new int [dim+1]; */
+
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY)
+ {
+ List_Element.push( ElNum );
+
+ //saving which nodes are boundary-nodes
+ for (int i=0; i<=dim; i++)
+ {
+ sv[i] = 1;
+ }
+ SMI_Get_elems(1,1,1,const_cast<DegreeOfFreedom*>(&ElNum),NULL,&nodeIndices,NULL,NULL);
+ SMI_Set_quantity_values(1,1,2,SMI_TYPE_INT,1,3,nodeIndices,sv);
+ }
+ //if the elemet isn't a boundary-element, but the interface cuts the FE in two nodes
+ else if (NumVertIntPoints == dim)
+ {
+ for (int i=0; i<=dim; i++)
+ {
+ if (elVertStatusVec[i] == ElementLevelSet::LEVEL_SET_BOUNDARY && elVertStatusVec[(i+1)%3] ==ElementLevelSet::LEVEL_SET_BOUNDARY)
+ {
+ value = 1;
+ SMI_Get_elems(1,1,1,const_cast<DegreeOfFreedom*>(&ElNum),NULL,&nodeIndices,NULL,NULL);
+ SMI_Set_quantity_values(1,1,2,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&nodeIndices[i]),&value);
+ SMI_Set_quantity_values(1,1,2,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&nodeIndices[(i+1)%3]),&value);
+
+ //saving which node of this element will be included into the second list
+ //part 1
+ for (int j=0; j<=dim; j++) //ACHTUNG: aus < ein <= gemacht
+ {
+ sv[j] = 0;
+ }
+ sv[(i+2)%3] = 1;
+
+ //include pair into the second list
+ vertStruct.ElNum = ElNum ;
+ vertStruct.VertNum = locInd[(i+2)%3];
+ List_El_Vert.push( vertStruct );
+
+ //saving which node of this element will be included into the second list
+ //part 2
+ SMI_Set_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,const_cast<DegreeOfFreedom*>(&ElNum),sv);
+
+ //counts how often a node is saved in the second list
+ if(count_in_list == 1)
+ {
+ SMI_Get_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct.VertNum),&value_q5);
+ value_q5 = value_q5 + 1;
+ SMI_Set_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct.VertNum),&value_q5);
+ }
+ }
+ }
+ }
+ else
+ {
+ //saving that none of the nodes of this element is saved in the second list
+ for (int i=0; i<=dim; i++)
+ {
+ sv[i] = 0;
+ }
+ SMI_Set_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,const_cast<DegreeOfFreedom*>(&ElNum),sv);
+ }
+ delete [] sv;
+ };
+
+//End function for creating two lists
+//
+
+//function for traversing the list "List_Element"
+//
+
+ void traverseListElement()
+ {
+ Vert_Struct vertStruct;
+
+ int Element;
+ int *neighbour=new int[dim+1];
+ int *oppVertices=new int[dim+1];
+ int *values=new int[dim*(dim+1)];
+ int *value = new int [dim+1];
+ int value_q2;
+ int value_q5;
+ int *nodeIndices;
+
+ while (List_Element.size() != 0)
+ {
+ Element = List_Element.front();
+ smiAdapter->getNeighbourInfo(Element, neighbour, oppVertices);
+
+ for ( int i=0; i<=dim; i++)
+ {
+ SMI_Get_elems(1,1,1,&neighbour[i],NULL,&nodeIndices, NULL,NULL);
+ SMI_Get_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,&neighbour[i],value);
+ SMI_Get_quantity_values(1,1,2,SMI_TYPE_INT,1,1,&nodeIndices[oppVertices[i]],&value_q2);
+ if (value[oppVertices[i]] == 0 && value_q2 == 0)
+ {
+ vertStruct.ElNum = neighbour[i] ;
+ vertStruct.VertNum = nodeIndices[oppVertices[i]] ;
+ List_El_Vert.push( vertStruct );
+ value[oppVertices[i]] = 1;
+
+ //counts how often a node is saved in the second list
+ if(count_in_list == 1)
+ {
+ SMI_Get_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct.VertNum),&value_q5);
+ value_q5 = value_q5 + 1;
+ SMI_Set_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct.VertNum),&value_q5);
+ }
+ }
+ SMI_Set_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,&neighbour[i],value);
+ }
+ List_Element.pop();
+ }
+ delete [] neighbour;
+ delete [] oppVertices;
+ delete [] values;
+ delete [] value;
+ };
+
+ //end of the function for traversing the list "List_Element"
+ //
+
+
+ /*
+ *gets the neighbour according to the node "VertNum_in" an two of its nodes
+ *returns "false" if the called neighbour exists and "true" if not
+ */
+
+ bool getNextNeighbour (const int ElNum_in, const int Vert_Up_in, const int VertNum_in, int &ElNum_out, int &VertNum_1_out, int &VertNum_2_out)
+ {
+ int *neighbour = new int[dim+1];
+ int *oppVertices = new int[dim+1];
+ int *nodeIndicesOfElem;
+ int locVertNum;
+
+ smiAdapter->getNeighbourInfo(ElNum_in, neighbour,oppVertices);
+
+ //which local node is the node "VertNum_in"?
+ SMI_Get_elems(1,1,1,const_cast<DegreeOfFreedom*>(&(ElNum_in)),NULL,&nodeIndicesOfElem, NULL, NULL);
+ for (int i=0; i<=dim; i++)
+ {
+ if(nodeIndicesOfElem[i] == VertNum_in)
+ {
+ locVertNum = i;
+ }
+ }
+
+ ElNum_out = neighbour [locVertNum];
+
+
+
+ SMI_Get_elems(1,1,1,const_cast<DegreeOfFreedom*>(&(ElNum_out)),NULL,&nodeIndicesOfElem, NULL, NULL);
+ VertNum_1_out = nodeIndicesOfElem [oppVertices[locVertNum]];
+ for (int i=0; i<=dim; i++)
+ {
+ if (nodeIndicesOfElem[i] != Vert_Up_in && nodeIndicesOfElem[i] != VertNum_1_out)
+ {
+ VertNum_2_out = nodeIndicesOfElem[i];
+ }
+ }
+
+ delete [] neighbour;
+ delete [] oppVertices;
+
+ if (ElNum_out == -1)
+ {
+ return true;
+ }
+
+ return false;
+
+ }
+
+
+ /*
+ * checking whetherthe element "ElNum_in" has to be included into the second list
+ *if yes it will be included
+ */
+
+ void includeIntoList (int ElNum_in, int VertNum_1_in, int VertNum_2_in)
+ {
+ Vert_Struct vertStruct_0;
+ Vert_Struct vertStruct_1;
+
+ int *nodeIndices;
+ int *valuesINT = new int[dim+1];
+ int *locVert = new int[dim];
+ int value_q2;
+ int value_q5;
+
+ if(ElNum_in != -1)
+ {
+
+ SMI_Get_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,const_cast<DegreeOfFreedom*>(&ElNum_in),valuesINT);
+
+ //which locla node is the node VerNum_1_in, which one is the node VertNum_2_in?
+ SMI_Get_elems (1, 1, 1, &ElNum_in, NULL, &nodeIndices, NULL, NULL);
+ for (int i=0; i<=dim; i++)
+ {
+ if(nodeIndices[i] == VertNum_1_in)
+ {
+ locVert[0] = i;
+ }
+ if(nodeIndices[i] == VertNum_2_in)
+ {
+ locVert[1] = i;
+ }
+ }
+
+ //if the pair of element and node isn't saved in the second List an if it isn't a boundary-node include into second list
+ SMI_Get_quantity_values(1,1,2,SMI_TYPE_INT,1,1,&VertNum_1_in, &value_q2);
+ if(valuesINT[locVert[0]] == 0 && value_q2 == 0)
+ {
+ vertStruct_0.ElNum = ElNum_in;
+ vertStruct_0.VertNum = VertNum_1_in;
+ List_El_Vert.push( vertStruct_0 );
+
+ valuesINT[locVert[0]] = 1;
+ SMI_Set_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,&ElNum_in,&valuesINT);
+
+ //counts how often a node is saved in the second list
+ if(count_in_list == 1)
+ {
+ SMI_Get_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct_0.VertNum),&value_q5);
+ value_q5 = value_q5 + 1;
+ SMI_Set_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct_0.VertNum),&value_q5);
+ }
+ }
+
+ //if the pair of element and node isn't saved in the second List an if it isn't a boundary-node include into second list
+ SMI_Get_quantity_values(1,1,2,SMI_TYPE_INT,1,1,&VertNum_2_in, &value_q2);
+ if(valuesINT[locVert[1]] == 0 && value_q2 == 0)
+ {
+ vertStruct_1.ElNum = ElNum_in;
+ vertStruct_1.VertNum = VertNum_2_in;
+ List_El_Vert.push( vertStruct_1 );
+
+ valuesINT[locVert[1]] = 1;
+ SMI_Set_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,&ElNum_in,&valuesINT);
+
+ //counts how often a node is saved in the second list
+ if(count_in_list == 1)
+ {
+ SMI_Get_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct_1.VertNum),&value_q5);
+ value_q5 = value_q5 + 1;
+ SMI_Set_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&vertStruct_1.VertNum),&value_q5);
+ }
+ }
+ }
+
+ delete [] valuesINT;
+ delete [] locVert;
+ }
+
+
+ /**
+ *function needed in the function "traversingListELVert"
+ **/
+ int getNext_node_l_r (int elem_l_r_in, int node_l_r_in, const int Vert)
+ {
+ int *nodeIndicesOfElem;
+
+ SMI_Get_elems(1,1,1,const_cast<DegreeOfFreedom*>(&elem_l_r_in),NULL,&nodeIndicesOfElem, NULL, NULL);
+ for (int i=0; i<=dim; i++)
+ {
+ if(nodeIndicesOfElem[i] != Vert && nodeIndicesOfElem[i] != node_l_r_in)
+ {
+ return nodeIndicesOfElem[i];
+ }
+ }
+ }
+
+ //Beginn function for traversing the list "List_El_Vert"
+ //
+
+ void traversingListELVert ( DOFVector<double> *boundVal_DOF )
+ {
+ Vert_Struct vertStruct_0;
+ Vert_Struct vertStruct_1;
+
+ Vert_Struct El_Vert;
+
+ int Vert;
+
+
+
+ //int * new_neighbour_l = new int [dim+1];
+ //int *new_neighbour_r = new int [dim+1];
+ // int *oppVertices_l = new int[dim+1];
+ //int *oppVertices_r = new int[dim+1];
+
+ int locVert;
+ int value_q3;
+ int value_q4;
+ int value_q5;
+ int counter = 0;
+ int *nodeIndices;
+ double *coords = new double [(dim+1)*dim];
+ int *valuesINT=new int[dim+1];
+ double update;
+ int value_q2;
+ int counter_list = 0;
+ FixVec<WorldVector<double> *, VERTEX> coordsOfNodes(dim);
+ for ( int i=0; i<=dim; i++)
+ {
+ coordsOfNodes[i] = new WorldVector<double>;
+ }
+ FixVec<double,VERTEX> boundVal(dim,NO_INIT);
+ WorldVector<double> helpVec;
+
+
+
+ while (List_El_Vert.size() != 0)
+ {
+ El_Vert = List_El_Vert.front();
+ Vert = El_Vert.VertNum;
+
+ SMI_Get_elems (1, 1, 1, &(El_Vert.ElNum), NULL, &(nodeIndices), NULL, NULL);
+ SMI_Get_nodes (1, 1, 3, dim, nodeIndices, coords);
+
+ counter = 0;
+ for (int i=0; i<=dim; i++)
+ {
+ if (nodeIndices[i] == El_Vert.VertNum)
+ {
+ (*(coordsOfNodes[dim]))[0]=coords[i*2];
+ (*(coordsOfNodes[dim]))[1]=coords[i*2+1];
+ boundVal[dim] = (*boundVal_DOF)[nodeIndices[i]];
+ locVert = i;
+ }
+ else
+ {
+ (*(coordsOfNodes[counter]))[0]=coords[i*2];
+ (*(coordsOfNodes[counter]))[1]=coords[i*2+1];
+ boundVal[counter] = (*boundVal_DOF)[nodeIndices[i]];
+ counter = counter+1;
+ }
+ }
+
+ update = elUpdate->calcElementUpdate( coordsOfNodes, boundVal );
+
+ if(count_updates == 1)
+ {
+ SMI_Get_quantity_values(1,1,3,SMI_TYPE_INT,1,1,&nodeIndices[locVert],&value_q3);
+ value_q3 = value_q3 + 1;
+ SMI_Set_quantity_values(1,1,3,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&nodeIndices[locVert]),&value_q3);
+ }
+
+ if (update < (*boundVal_DOF)[El_Vert.VertNum] && ((*boundVal_DOF)[El_Vert.VertNum]-update) > tol)
+ {
+ (*boundVal_DOF)[El_Vert.VertNum] = update;
+
+ if(count_updates == 1)
+ {
+ SMI_Get_quantity_values(1,1,4,SMI_TYPE_INT,1,1,&nodeIndices[locVert],&value_q4);
+ value_q4 = value_q4 + 1;
+ SMI_Set_quantity_values(1,1,4,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(&nodeIndices[locVert]),&value_q4);
+ }
+ search_and_include_comb(El_Vert.ElNum, Vert, nodeIndices);
+ }
+
+ SMI_Get_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,&El_Vert.ElNum, valuesINT);
+ List_El_Vert.pop();
+ valuesINT[locVert] = 0;
+ SMI_Set_quantity_values(1,1,1,SMI_TYPE_INT,dim+1,1,&El_Vert.ElNum, valuesINT);
+
+ //counts how often a node is saved in the second list
+ if(count_in_list == 1)
+ {
+ SMI_Get_quantity_values(1,1,5,SMI_TYPE_INT,1,1,const_cast<DegreeOfFreedom*>(& El_Vert.VertNum),&value_q5);
+ value_q5 = value_q5 - 1;
+ SMI_Set_quantity_values(1,1,5,SMI_TYPE_INT,dim+1,1,const_cast<DegreeOfFreedom*>(& El_Vert.VertNum),&value_q5);
+ }
+
+ if(count_in_list == 1)
+ {
+ counter_list++;
+ if(counter_list % print_in_list == 0)
+ {
+ print_quantity_5(counter_list);
+ }
+ }
+ }
+
+// delete [] new_neighbour_l;
+// delete [] new_neighbour_r;
+// delete [] oppVertices_l;
+// delete [] oppVertices_r;
+ delete [] coords;
+ delete [] valuesINT;
+ for ( int i=0; i<Global::getGeo(VERTEX,dim); i++)
+ {
+ delete coordsOfNodes[i];
+ }
+ };
+
+ //End of the function for traversing the list "List_El_Vert"
+ //
+
+ void search_and_include_comb(int ElNum, int Vert, int *nodeIndices)
+ {
+ bool stop_l = false;
+ bool stop_r = false;
+ int neighbour_l;
+ int neighbour_r;
+ int Vert_1_l;
+ int Vert_2_l;
+ int Vert_1_r;
+ int Vert_2_r;
+ int locVert;
+ int counter;
+ int elem_l;
+ int elem_r;
+ int node_l;
+ int node_r;
+
+ //here the process of including elements/nodes into the second list
+ stop_l = false;
+ stop_r = false;
+
+ elem_l = ElNum;
+ if (elem_l == -1)
+ {
+ stop_l = true;
+ }
+ elem_r = ElNum;
+ if (elem_r == -1)
+ {
+ stop_r = true;
+ }
+
+ counter = 0;
+ for (int i=0; i<=dim; i++)
+ {
+ if(nodeIndices[i] != Vert && counter == 1)
+ {
+ node_r = nodeIndices[i];
+ counter = counter + 1;
+ }
+ if(nodeIndices[i] != Vert && counter == 0)
+ {
+ node_l = nodeIndices[i];
+ counter = counter + 1;
+ }
+
+ if(nodeIndices[i] == Vert)
+ {
+ locVert = i;
+ }
+ }
+
+ while (!stop_l || !stop_r)
+ {
+ //get next neighbour
+ if (!stop_l)
+ {
+ stop_l = getNextNeighbour (elem_l, Vert, node_l, neighbour_l, Vert_1_l, Vert_2_l);
+ }
+ if(!stop_r)
+ {
+ stop_r = getNextNeighbour (elem_r, Vert, node_r, neighbour_r, Vert_1_r, Vert_2_r);
+ }
+
+ if(neighbour_l == elem_r)
+ {
+ break;
+ }
+
+ //indclude into the second list (only if possible)
+ if (!stop_l)
+ {
+ includeIntoList (neighbour_l, Vert_1_l, Vert_2_l);
+ if(neighbour_l == neighbour_r)
+ {
+ break;
+ }
+ //"elem_l", "node_l", "elem_r", "node_r" have to be set on the next elements
+ node_l = getNext_node_l_r(elem_l, node_l,Vert);
+ elem_l = neighbour_l;
+ }
+ if (!stop_r && neighbour_l != neighbour_r)
+ {
+ includeIntoList (neighbour_r, Vert_1_r, Vert_2_r);
+ //"elem_l", "node_l", "elem_r", "node_r" have to be set on the next elements
+ node_r = getNext_node_l_r(elem_r, node_r, Vert);
+ elem_r = neighbour_r;
+ }
+ }
+ }
+
+ //==================================
+ void print_quantity_5 (int cntr)
+ {
+ int numNodes;
+ int*nodeIndices;
+ double *coords = new double[dim];
+ SMI_Get_all_nodes(1,1,const_cast<DegreeOfFreedom*>(&numNodes),&nodeIndices);
+ int *values_q5 = new int [numNodes];
+
+ std::string q5_OutFile_0;
+ std::string q5_OutFile_1;
+ std::string q5_OutFile_2;
+ std::string q5_OutFile_3;
+ std::string q5_OutFile_4;
+
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_0", &q5_OutFile_0);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_1", &q5_OutFile_1);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_2", &q5_OutFile_2);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_3", &q5_OutFile_3);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_4", &q5_OutFile_4);
+
+
+ char cntrStr[20];
+ sprintf(cntrStr, "%d", cntr);
+ q5_OutFile_0 += cntrStr;
+ q5_OutFile_1 += cntrStr;
+ q5_OutFile_2 += cntrStr;
+ q5_OutFile_3 += cntrStr;
+ q5_OutFile_4 += cntrStr;
+
+ cout << "count saving Ausgabe-Datei_0: " << q5_OutFile_0.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_1: " << q5_OutFile_1.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_2: " << q5_OutFile_2.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_3: " << q5_OutFile_3.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_4: " << q5_OutFile_4.c_str() << "\n\n";
+
+ ofstream out_0 (q5_OutFile_0.c_str());
+ ofstream out_1 (q5_OutFile_1.c_str());
+ ofstream out_2 (q5_OutFile_2.c_str());
+ ofstream out_3 (q5_OutFile_3.c_str());
+ ofstream out_4 (q5_OutFile_4.c_str());
+
+ SMI_Get_quantity_values(1,1,5,SMI_TYPE_INT,1,numNodes,nodeIndices,values_q5);
+
+ for (int i=0; i<numNodes; i++)
+ {
+
+ SMI_Get_nodes (1,1,1,dim,const_cast<DegreeOfFreedom*>(&nodeIndices[i]),coords);
+
+ if(values_q5[i] == 0)
+ {
+ out_0<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] == 1)
+ {
+ out_1<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] == 2)
+ {
+ out_2<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] == 3)
+ {
+ out_3<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] >= 4)
+ {
+ out_4<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+
+
+ }
+
+ out_0.close();
+ out_1.close();
+ out_2.close();
+ out_3.close();
+ out_4.close();
+
+ delete [] coords;
+ delete [] values_q5;
+ }
+ //====================================
+
+ protected:
+ /**
+ * Tolerance for Hopf-Lax update iteration loop.
+ */
+ double tol;
+
+ /**
+ *is needed for transfering the mesh to SMI
+ */
+ SMIAdapter *smiAdapter;
+
+ /**
+ * in this list boundary-elements are saved;
+ * is needed for creating the list "List_El_Vert"
+ */
+ queue<int> List_Element;
+
+ /**
+ *in this list structs filled with element-number and vertex-number are saved;
+ *is needed for traversing the mesh efficiently
+ */
+ queue<Vert_Struct> List_El_Vert;
+
+ /**
+ *0 ->do not count updates
+ *1 ->count updates
+ **/
+ int count_updates;
+
+ /**
+ *1 -> it will be count how often a node is saved in the second list
+ *0 ->it will not be count
+ **/
+ int count_in_list;
+
+ int print_in_list;
+};
+
+#endif // HL_SIGNEDDISTBORNEMANN
diff --git a/AMDiS/Reinit/src/HL_SignedDistLevels.cc b/AMDiS/Reinit/src/HL_SignedDistLevels.cc
new file mode 100644
index 00000000..4400e22d
--- /dev/null
+++ b/AMDiS/Reinit/src/HL_SignedDistLevels.cc
@@ -0,0 +1,1360 @@
+#include "HL_SignedDistLevels.h"
+
+#include "VelocityExtFromVelocityField.h"
+
+void
+HL_SignedDistLevels::initializeBoundary()
+{
+ FUNCNAME("HL_SignedDistLevels::initializeBoundary()");
+
+ TraverseStack stack;
+ FixVec<double, VERTEX> distVec(dim, NO_INIT);
+ int elStatus;
+ const int *elVertStatusVec;
+ int NumVertIntPoints;
+ int ElNum;
+ clock_t time1;
+ clock_t time2;
+ double TIME;
+
+ //=== transvering Mesh to SMI and add quantities
+ time1 = clock();
+ Mesh_to_SMI_and_quantity();
+ time2 = clock();
+ TIME = TIME_USED(time1, time2);
+
+ cout <<"Zeit zum Transformieren nach SMI: "<<TIME<<" sec.\n";
+
+ // ===== All non-boundary vertices are initialized with "infinity". =====
+ sD_DOF->set(inftyValue);
+
+ // ===== Traverse mesh and initialize boundary elements. =====
+ const int nBasFcts = feSpace->getBasisFcts()->getNumber();
+ DegreeOfFreedom *locInd = GET_MEMORY(DegreeOfFreedom, nBasFcts);
+
+ ElInfo *elInfo;
+ if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
+ elInfo = stack.traverseFirst(feSpace->getMesh(),
+ -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_GRD_LAMBDA);
+ }
+ else {
+ elInfo = stack.traverseFirst(feSpace->getMesh(),
+ -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_COORDS);
+ }
+ while(elInfo) {
+
+ // Set elInfo in case velocity extension from velocity field is used.
+ if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
+ ((VelocityExtFromVelocityField *)velExt)->setElInfo(elInfo);
+ }
+
+ // Get local indices of vertices.
+ feSpace->getBasisFcts()->getLocalIndices(
+ const_cast<Element *>(elInfo->getElement()),
+ const_cast<DOFAdmin *>(feSpace->getAdmin()),
+ locInd);
+
+ // Get element status.
+ elStatus = elLS->createElementLevelSet(elInfo);
+
+ //Get some information for creating the first list
+
+ elVertStatusVec = elLS->getElVertStatusVec();
+
+ //Beginn creating the first list
+
+ NumVertIntPoints = elLS->getNumVertIntPoints();
+ ElNum = elInfo->getElement()->getIndex();
+
+ createFirstList( elStatus, elVertStatusVec, ElNum, NumVertIntPoints);
+
+ //end creating the first list
+
+ // Is element cut by the interface ?
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+
+ // Reset element distance vector.
+ for (int i=0; i<=dim; ++i) {
+ distVec[i] = inftyValue;
+ }
+
+ // Mark all vertices as boundary vertices.
+ for (int i=0; i<=dim; ++i) {
+ (*bound_DOF)[locInd[i]] = 1.0;
+ }
+
+ // Calculate distance for all vertices.
+ if (bndElDist->calcDistOnBoundaryElement(elInfo, distVec) !=
+ ElementLevelSet::LEVEL_SET_BOUNDARY) {
+ ERROR_EXIT("error in distance calculation !\n");
+ }
+ else {
+
+ // If distance is smaller, correct to new distance.
+ for (int i=0; i<=dim; ++i) {
+ if ((*sD_DOF)[locInd[i]] > distVec[i]) {
+ (*sD_DOF)[locInd[i]] = distVec[i];
+ //If distance is corrected, calculate new velocity.
+ if(velExt != NULL)
+ {
+ velExt->calcVelocityBoundary(locInd, i);
+ }
+ }
+ }
+ }
+ } // end of: elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY
+
+ else if (elLS->getNumVertIntPoints() != 0){
+// // Interface cuts element only in vertices.
+// elVertStatusVec = elLS->getElVertStatusVec();
+
+// for (int i=0; i<=dim; ++i) {
+// if (elVertStatusVec[i] == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+// (*bound_DOF)[locInd[i]] = 1.0;
+// (*sD_DOF)[locInd[i]] = 0.0;
+// }
+// }
+
+ ERROR_EXIT("intersection point is element vertex. should never happen !\n");
+ }
+
+ elInfo = stack.traverseNext(elInfo);
+ } // end of: mesh traverse
+
+ FREE_MEMORY(locInd, DegreeOfFreedom, nBasFcts);
+
+ return;
+}
+
+void
+HL_SignedDistLevels::Mesh_to_SMI_and_quantity()
+{
+ if(!smiAdapter)
+ {
+ smiAdapter = NEW SMIAdapter(1, 1,
+ const_cast<FiniteElemSpace*>(feSpace),
+ 1, -1);
+ }
+
+ // cout << "\n\n\tSMI-Adapter angelegt !\n\n";
+
+ //====== transfer Mesh to SMI ======================
+
+ smiAdapter->addNeighbourInfo();
+
+ // cout << "\n\n\tNachbarschafts-Infos in SMI ergaenzt !\n\n";
+
+ smiAdapter-> transferMeshToSMI();
+
+ // cout << "\n\n\tGitter nach SMI geschrieben !\n\n";
+
+ int smiError = SMI_Begin_write_transaction(1,1);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Begin_write_transaction() failed with error %d\n", smiError);
+
+ int nul = 0;
+ int *null = new int [dim+1];
+ for (int i=0; i<=dim; i++)
+ {
+ null[i] = nul;
+ }
+
+ //which pair of element and node is saved in list
+ smiError =
+ SMI_Add_quantity(1, 1, 1, SMI_LOC_ELEM, SMI_TYPE_INT, dim+1, null);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Add_quantity() failed with error %d\n", smiError);
+
+ //which node is a boundary-node
+ smiError = SMI_Add_quantity(1, 1, 2, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Add_quantity() failed with error %d\n", smiError);
+
+ //saves the number of tried updates
+ smiError = SMI_Add_quantity(1, 1, 3, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Add_quantity() failed with error %d\n", smiError);
+
+ //saves the number of updates
+ smiError = SMI_Add_quantity(1, 1, 4, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Add_quantity() failed with error %d\n", smiError);
+
+ //how often a node is saved in the second list
+ smiError = SMI_Add_quantity(1, 1, 5, SMI_LOC_NODE, SMI_TYPE_INT, 1, &nul);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Add_quantity() failed with error %d\n", smiError);
+
+ //saves the level of an element
+ smiError = SMI_Add_quantity(1, 1, 6, SMI_LOC_ELEM, SMI_TYPE_INT, 1,
+ &inftyValue);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Add_quantity() failed with error %d\n", smiError);
+
+ // cout << "\n\n\tQuantities in SMI ergaenzt !\n\n";
+
+ delete [] null;
+}
+
+void
+HL_SignedDistLevels::createFirstList(const int elStatus,
+ const int *elVertStatusVec,
+ int ElNum,
+ const int NumVertIntPoints)
+{
+ int *sv = new int[dim+1];
+ int *nodeIndices;
+ int nul = 0;
+ int smiError;
+
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY)
+ {
+ List_Element.push(ElNum);
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, 1, &ElNum, &nul);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n", smiError);
+
+ //saving which nodes are boundary-nodes
+ for (int i=0; i<=dim; i++)
+ {
+ sv[i] = 1;
+ }
+
+ smiError =
+ SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&ElNum),
+ NULL, &nodeIndices, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 3, nodeIndices, sv);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n", smiError);
+ }
+ //if the elemet isn't a boundary-element, but the interface cuts the FE in two nodes
+// else if (NumVertIntPoints == dim)
+// {
+// List_Element.push( ElNum );
+
+// smiError =
+// SMI_Set_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, 1, &ElNum, &nul);
+// TEST_EXIT(smiError == SMI_OK)
+// ("SMI_Set_quantity_values() failed with error %d\n", smiError);
+
+// //saving which nodes are boundary-nodes
+// for (int i=0; i<=dim; i++)
+// {
+// if(elVertStatusVec[i] == ElementLevelSet::LEVEL_SET_BOUNDARY)
+// {
+// sv[i] = 1;
+// }
+// else
+// {
+// sv[i] = 0;
+// }
+
+// smiError =
+// SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&ElNum),
+// NULL, &nodeIndices, NULL, NULL);
+// TEST_EXIT(smiError == SMI_OK)
+// ("SMI_Get_elems() failed with error %d\n", smiError);
+
+// smiError =
+// SMI_Set_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 3,
+// nodeIndices, sv);
+// TEST_EXIT(smiError == SMI_OK)
+// ("SMI_Set_quantity_values() failed with error %d\n", smiError);
+// }
+// }
+
+ delete [] sv;
+}
+
+void
+HL_SignedDistLevels::HL_updateIteration()
+{
+ int numNodes;
+ int *nodeIndices;
+ int max_q3 = 0;
+ int max_q4 = 0;
+ int sum_q3 = 0;
+ int sum_q4 = 0;
+ clock_t time1;
+ clock_t time2;
+ double TIME;
+ int control;
+
+ int smiError =
+ SMI_Get_all_nodes(1, 1, const_cast<DegreeOfFreedom*>(&numNodes),
+ &nodeIndices);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_all_nodes() failed with error %d\n", smiError);
+
+ int *values_q3 = new int [numNodes];
+ int *values_q4 = new int [numNodes];
+
+ GET_PARAMETER(0,"SignedDist->count updates", "%d", &count_updates);
+
+ time1 = clock();
+ control = traverseListElement();
+ time2 = clock();
+ TIME = TIME_USED(time1, time2);
+
+ cout<<"Zeit zum durchlaufen der ersten (Elementen-) Liste: "<<TIME<<" sec.\n\n";
+ if(control == 1)
+ {
+ createLevels ();
+ }
+
+ /////////////////////////////////
+
+
+ //print chosen level in a file
+ if (print_level == 1)
+ {
+ DOFVector<double> *level_DOF = NEW DOFVector<double>(sD_DOF->getFESpace(),
+ "level_DOF");
+
+ int numElems;
+ int *elemIndices;
+ int *NodeIndices;
+
+ smiError = SMI_Get_all_elems(1, 1, &numElems, &elemIndices);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_all_elems() failed with error %d\n", smiError);
+
+ int *value_q6 = new int [numElems];
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, numElems,
+ elemIndices, value_q6);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ for(int i=0; i<numElems; i++)
+ {
+ smiError =
+ SMI_Get_elems(1, 1, 1, &elemIndices[i], NULL, &NodeIndices,
+ NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ for(int j=0; j<=dim; j++)
+ {
+ (*level_DOF)[NodeIndices[j]] = -1;
+ }
+ }
+ for(int i=0; i<numElems; i++)
+ {
+ smiError =
+ SMI_Get_elems(1, 1, 1, &elemIndices[i], NULL, &NodeIndices,
+ NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ if (value_q6[i] == chosen_level)
+ {
+ for(int j=0; j<dim+1; j++)
+ {
+ (*level_DOF)[NodeIndices[j]] = 0;
+ }
+ }
+ }
+
+ FileWriter *levelFileWriter =
+ NEW FileWriter("level->output",
+ level_DOF->getFESpace()->getMesh(),
+ level_DOF);
+ levelFileWriter->writeFiles (adaptInfo, false);
+
+ DELETE levelFileWriter;
+ }
+
+ //end of the part for printing chosen levels
+
+
+ time1 = clock();
+ traversingListLevel_ ( sD_DOF );
+ time2 = clock();
+ TIME = TIME_USED(time1,time2);
+
+ cout<<"Zeit zum Durchlaufen der zweiten Liste (Liste von Listen): "<<TIME<<" sec.\n\n";
+
+ std::string smiOutFile;
+ GET_PARAMETER(0, "SignedDist->count updates->output->filename", &smiOutFile);
+ cout << "count updates Ausgabe-Datei: " << smiOutFile.c_str() << "\n\n";
+
+ ofstream out (smiOutFile.c_str());
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 3, SMI_TYPE_INT, 1, numNodes, nodeIndices,
+ values_q3);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 4, SMI_TYPE_INT, 1, numNodes, nodeIndices,
+ values_q4);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ for (int i=0; i<numNodes; i++)
+ {
+ out<<nodeIndices[i]<<" "<<values_q3[i]<<" "<<values_q4[i]<<"\n";
+
+ if(max_q3 < values_q3[i])
+ {
+ max_q3 = values_q3[i];
+ }
+ if(max_q4 < values_q4[i])
+ {
+ max_q4 = values_q4[i];
+ }
+ sum_q3 = sum_q3 + values_q3[i];
+ sum_q4 = sum_q4 + values_q4[i];
+ }
+
+ out<<"\n\n maximale Anzahl an versuchten Updates auf einem Knoten: "<<max_q3<<"\n maximale Anzahl an durchgefuehrten Updates auf einem Knoten: "<<max_q4<<"\n";
+ out<<"\n Summe aller versuchten Updates: "<<sum_q3<<"\n Summe aller durchgefuehrten updates: "<<sum_q4<<"\n";
+
+ out<<"Anzahl an Knoten: "<<numNodes<<"\n";
+
+ out.close();
+
+ smiError = SMI_End_write_transaction(1, 1);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_End_write_transaction() failed with error %d\n", smiError);
+
+ return;
+}
+
+int
+HL_SignedDistLevels::traverseListElement()
+{
+ int Element;
+
+ bool neighboursExist;
+ bool neighbourInNewListSet = false;
+
+ while (List_Element.size() != 0)
+ {
+ Element = List_Element.front();
+
+ neighboursExist = collectNeighbours_setLevels (Element, 0, &neighbourInNewListSet);
+
+ List_Element.pop();
+ }
+
+ //within the function "collectNeighbours_setLevels" the new pairs are saved in the list "level"
+ //this list will be one entry of the list "Level_"
+ if (neighbourInNewListSet)
+ {
+ Level_.push_back (Level);
+ }
+
+ while (Level.size() != 0)
+ {
+ Level.pop();
+ }
+ if (!neighbourInNewListSet)
+ {
+ return 0;
+ }
+ else
+ {
+ return 1;
+ }
+
+}
+
+void
+HL_SignedDistLevels::createLevels ()
+{
+ bool neighbourExists = true;
+ bool stop = false;
+ Vert_Struct vertStruct;
+ int i = 0;
+ bool elementInNextListSet = true;
+ int value_q2;
+ int smiError;
+
+ while(elementInNextListSet)
+ {
+ stop = false;
+ neighbourExists = false;
+ elementInNextListSet = false;
+
+ while(Level_[i].size() != 0)
+ {
+ vertStruct = Level_[i].front();
+
+ neighbourExists = collectNeighbours_setLevels (vertStruct.ElNum, i+1, &elementInNextListSet);
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 1,
+ &vertStruct.VertNum, &value_q2);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ if(value_q2 == 0)
+ {
+ helpLevel.push (vertStruct);
+ }
+ Level_[i].pop();
+ }
+
+ Level_[i] = helpLevel;
+ while (helpLevel.size() != 0)
+ {
+ helpLevel.pop();
+ }
+
+ if(elementInNextListSet)
+ {
+ Level_.push_back (Level);
+ }
+
+ while (Level.size() != 0)
+ {
+ Level.pop();
+ }
+ i++;
+ }
+}
+
+bool
+HL_SignedDistLevels::collectNeighbours_setLevels (const int Element,
+ const int currentIndex,
+ bool *elementInNewListSet)
+{
+ Vert_Struct vertStruct;
+
+ int *neighbour=new int[dim+1];
+ int *oppVertices=new int[dim+1];
+ int *value = new int [dim+1];
+ int saved_level;
+ int value_q5;
+ int *nodeIndices;
+ int *nodeIndicesOfElem;
+ int smiError;
+
+ bool neighbourExists = false;
+
+ smiAdapter->getNeighbourInfo(Element, neighbour, oppVertices);
+
+ for ( int i=0; i<=dim; i++)
+ {
+ //if this neighbour exists it will be saved in a list according to it's level
+ if(neighbour[i] != -1)
+ {
+ neighbourExists = true;
+
+ //calling old level and saving new level if necessary
+ smiError =
+ SMI_Get_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, 1,
+ &neighbour[i], &saved_level);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ if(saved_level > currentIndex)
+ {
+ saved_level = currentIndex+1;
+
+ smiError =
+ SMI_Set_quantity_values (1, 1, 6, SMI_TYPE_INT, 1, 1,
+ &neighbour[i], &saved_level);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+
+ smiError =
+ SMI_Get_elems(1, 1, 1, &neighbour[i], NULL, &nodeIndices,
+ NULL,NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n",
+ smiError);
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 1, SMI_TYPE_INT, dim+1, 1,
+ &neighbour[i], value);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ if (value[oppVertices[i]] == 0)
+ {
+ vertStruct.ElNum = neighbour[i] ;
+ vertStruct.VertNum = oppVertices[i];
+ Level.push( vertStruct );
+ value[oppVertices[i]] = 1;
+ (*elementInNewListSet) = true;
+
+ //counts how often a node is saved in the list of lists
+ if(count_in_list == 1)
+ {
+ smiError =
+ SMI_Get_elems(1, 1, 1,
+ const_cast<DegreeOfFreedom*>(&(vertStruct.ElNum)),
+ NULL, &nodeIndicesOfElem, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndicesOfElem[vertStruct.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ value_q5 = value_q5 + 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndicesOfElem[vertStruct.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+ }
+ }
+ //save which pairs of elements and nodes is saved in a list
+ smiError =
+ SMI_Set_quantity_values(1, 1, 1, SMI_TYPE_INT, dim+1, 1,
+// &neighbour[i], value);
+ neighbour+i, value);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+ }
+ }
+ }
+
+ delete [] neighbour;
+ delete [] oppVertices;
+ delete [] value;
+
+ return neighbourExists;
+
+
+}
+
+void
+HL_SignedDistLevels::traversingListLevel_ ( DOFVector<double> *boundVal_DOF )
+{
+ bool Continue = true;
+
+ Vert_Struct El_Vert;
+
+ int Vert;
+
+ int locVert;
+ int value_q3;
+ int value_q4;
+ int value_q5;
+ int counter = 0;
+ int *nodeIndices;
+ int *nodeIndicesOfElem;
+ double *coords = new double [(dim+1)*dim];
+ int *valuesINT=new int[dim+1];
+ double update;
+ int counter_list = 0;
+ int counter_lists_of_lists = 0;
+ FixVec<WorldVector<double> *, VERTEX> coordsOfNodes(dim);
+ for ( int i=0; i<=dim; i++)
+ {
+ coordsOfNodes[i] = new WorldVector<double>;
+ }
+ FixVec<double,VERTEX> boundVal(dim,NO_INIT);
+ WorldVector<double> helpVec;
+ int smiError;
+
+ //for repeating the traverse of the list "Level_"
+ while (Continue == true)
+ {
+ Continue = false;
+
+ //for traversing of the list "Level_"
+ for (int j=0; j < (signed int) Level_.size();j++)
+ {
+ //for traversing the entries of the List "Level_"
+
+ while (Level_[j].size() != 0)
+ {
+ El_Vert = Level_[j].front();
+ Vert = El_Vert.VertNum;
+
+ smiError =
+ SMI_Get_elems (1, 1, 1, &(El_Vert.ElNum), NULL,
+ &(nodeIndices), NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ smiError =
+ SMI_Get_nodes (1, 1, 3, dim, nodeIndices, coords);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_nodes() failed with error %d\n", smiError);
+
+ //sorting of the nodes for calculating the update
+ counter = 0;
+ for (int i=0; i<=dim; i++)
+ {
+ if (i == El_Vert.VertNum)
+ {
+ (*(coordsOfNodes[dim]))[0]=coords[i*2];
+ (*(coordsOfNodes[dim]))[1]=coords[i*2+1];
+ boundVal[dim] = (*boundVal_DOF)[nodeIndices[i]];
+ locVert = i;
+ }
+ else
+ {
+ (*(coordsOfNodes[counter]))[0]=coords[i*2];
+ (*(coordsOfNodes[counter]))[1]=coords[i*2+1];
+ boundVal[counter] = (*boundVal_DOF)[nodeIndices[i]];
+ counter = counter+1;
+ }
+ }
+
+ //save permutation of vertexes for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setPermutation(locVert, 0);
+ }
+ update = elUpdate->calcElementUpdate( coordsOfNodes, boundVal );
+
+ //for counting tried updates
+ if(count_updates == 1)
+ {
+ smiError =
+ SMI_Get_quantity_values(1, 1, 3, SMI_TYPE_INT, 1, 1,
+ &nodeIndices[locVert], &value_q3);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ value_q3 = value_q3 + 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 3, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(&nodeIndices[locVert]),
+ &value_q3);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+ }
+
+ //checking whether the calculated update will be set as new distance
+ if (update < (*boundVal_DOF)[nodeIndices[locVert]] && ((*boundVal_DOF)[nodeIndices[locVert]]-update) > tol)
+ {
+ (*boundVal_DOF)[nodeIndices[locVert]] = update;
+ //If distance is corrected, calculate new velocity.
+ if(velExt != NULL)
+ {
+ velExt->calcVelocity(nodeIndices, locVert);
+ }
+
+ //for counting successful updates
+ if(count_updates == 1)
+ {
+ smiError =
+ SMI_Get_quantity_values(1, 1, 4, SMI_TYPE_INT, 1, 1,
+ &nodeIndices[locVert],
+ &value_q4);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ value_q4 = value_q4 + 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 4, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(&nodeIndices[locVert]),
+ &value_q4);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+ }
+ //??
+ Continue = search_and_include_comb(El_Vert.ElNum,
+ Vert, nodeIndices, j);
+ }
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 1, SMI_TYPE_INT, dim+1, 1,
+ &El_Vert.ElNum, valuesINT);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ Level_[j].pop();
+ valuesINT[locVert] = 0;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 1, SMI_TYPE_INT, dim+1, 1,
+ &El_Vert.ElNum, valuesINT);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+
+ //counts how often a node is saved in the list of lists
+ if(count_in_list == 1)
+ {
+ smiError =
+ SMI_Get_elems(1, 1, 1,
+ const_cast<DegreeOfFreedom*>(&(El_Vert.ElNum)),
+ NULL, &nodeIndicesOfElem, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndicesOfElem[El_Vert.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ value_q5 = value_q5 - 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndicesOfElem[El_Vert.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+ }
+
+ //for printing which node is how often saved in a list after a special number of tried updates
+ if(count_in_list == 1)
+ {
+ counter_list++;
+ if(counter_list % print_in_list == 0)
+ {
+ print_quantity_5(counter_list);
+ }
+ }
+ }
+ }
+ //for printing which node is how often saved in a list at the end of the ..th traverse of the list "Level_"
+ if(count_in_list == 1)
+ {
+ counter_lists_of_lists++;
+ if(counter_lists_of_lists == print_in_list_2)
+ {
+ print_quantity_5(counter_lists_of_lists);
+ }
+ }
+
+ }
+
+ delete [] coords;
+ delete [] valuesINT;
+ for ( int i=0; i<Global::getGeo(VERTEX,dim); i++)
+ {
+ delete coordsOfNodes[i];
+ }
+}
+
+bool
+HL_SignedDistLevels::search_and_include_comb(int ElNum,
+ int Vert,
+ int *nodeIndices,
+ const int Index)
+{
+ bool stop_l = false;
+ bool stop_r = false;
+ bool Continue = false;
+ bool continue_out = false;
+ int neighbour_l;
+ int neighbour_r;
+ int Vert_1_l;
+ int Vert_2_l;
+ int Vert_1_r;
+ int Vert_2_r;
+ int counter;
+ int elem_l;
+ int elem_r;
+ int node_l;
+ int node_r;
+ int VertGlobal;
+ int *nodeIndicesOfElem;
+
+ int smiError =
+ SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&(ElNum)), NULL,
+ &nodeIndicesOfElem, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ VertGlobal = nodeIndicesOfElem[Vert];
+
+ //here the process of including elements/nodes into the second list
+ stop_l = false;
+ stop_r = false;
+
+ elem_l = ElNum;
+ if (elem_l == -1)
+ {
+ stop_l = true;
+ }
+ elem_r = ElNum;
+ if (elem_r == -1)
+ {
+ stop_r = true;
+ }
+
+ counter = 0;
+ for (int i=0; i<=dim; i++)
+ {
+ if(i != Vert && counter == 1)
+ {
+ node_r = i;
+ counter = counter + 1;
+ }
+ if(i != Vert && counter == 0)
+ {
+ node_l = i;
+ counter = counter + 1;
+ }
+
+ }
+
+ while (!stop_l || !stop_r)
+ {
+ //get next neighbour
+ if (!stop_l)
+ {
+ stop_l = getNextNeighbour (elem_l, VertGlobal, node_l, neighbour_l, Vert_1_l, Vert_2_l);
+ }
+ if(!stop_r)
+ {
+ stop_r = getNextNeighbour (elem_r, VertGlobal, node_r, neighbour_r, Vert_1_r, Vert_2_r);
+ }
+
+ if(neighbour_l == elem_r)
+ {
+ break;
+ }
+
+ //indclude into the second list (only if possible)
+ if (!stop_l)
+ {
+ Continue = includeIntoList (neighbour_l, Vert_1_l, Vert_2_l, Index);
+ if(Continue == true)
+ {
+ continue_out = Continue;
+ }
+ if(neighbour_l == neighbour_r)
+ {
+ break;
+ }
+ //"elem_l", "node_l", "elem_r", "node_r" have to be set on the next elements
+ node_l = getNext_node_l_r(elem_l,neighbour_l, node_l,VertGlobal);
+ elem_l = neighbour_l;
+
+ }
+ if (!stop_r && neighbour_l != neighbour_r)
+ {
+ Continue = includeIntoList (neighbour_r, Vert_1_r, Vert_2_r, Index);
+ if(Continue == true)
+ {
+ continue_out = Continue;
+ }
+ //"elem_l", "node_l", "elem_r", "node_r" have to be set on the next elements
+ node_r = getNext_node_l_r(elem_r,neighbour_r, node_r, VertGlobal);
+ elem_r = neighbour_r;
+ }
+ }
+ return continue_out;
+}
+
+bool
+HL_SignedDistLevels::includeIntoList (int ElNum_in,
+ int VertNum_1_in,
+ int VertNum_2_in,
+ const int Index)
+{
+ Vert_Struct vertStruct_0;
+ Vert_Struct vertStruct_1;
+
+ int *nodeIndices;
+ int *valuesINT = new int[dim+1];
+ int *locVert = new int[dim];
+ int value_q2;
+ int value_q5;
+ int value_q6;
+ int smiError;
+
+ bool Continue = false;
+
+ if(ElNum_in != -1)
+ {
+ smiError =
+ SMI_Get_elems(1, 1, 1, &ElNum_in, NULL, &nodeIndices, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 1, SMI_TYPE_INT, dim+1, 1,
+ const_cast<DegreeOfFreedom*>(&ElNum_in),
+ valuesINT);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ //if the pair of element and node isn't saved in any List and if it isn't a boundary-node include into adapting list
+ smiError = SMI_Get_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 1,
+ &nodeIndices[VertNum_1_in],
+ &value_q2);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ if(valuesINT[VertNum_1_in] == 0 && value_q2 == 0)
+ {
+ vertStruct_0.ElNum = ElNum_in;
+ vertStruct_0.VertNum = VertNum_1_in;
+
+ smiError =
+ SMI_Get_quantity_values (1, 1, 6, SMI_TYPE_INT, 1, 1, &ElNum_in,
+ &value_q6);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ Level_[value_q6-1].push( vertStruct_0 );
+
+ //counts how often a node is saved in the list of lists
+ if(count_in_list == 1)
+ {
+ smiError =
+ SMI_Get_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndices[vertStruct_0.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ value_q5 = value_q5 + 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndices[vertStruct_0.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+ }
+
+ if(value_q6-1 < Index)
+ {
+ Continue = true;
+ }
+
+ valuesINT[VertNum_1_in] = 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 1, SMI_TYPE_INT, dim+1, 1,
+// &ElNum_in, &valuesINT);
+ &ElNum_in, valuesINT);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n", smiError);
+ }
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 2, SMI_TYPE_INT, 1, 1,
+ &nodeIndices[VertNum_2_in], &value_q2);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ if(valuesINT[VertNum_2_in] == 0 && value_q2 == 0)
+ {
+ vertStruct_1.ElNum = ElNum_in;
+ vertStruct_1.VertNum = VertNum_2_in;
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, 1, &ElNum_in,
+ &value_q6);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ Level_[value_q6-1].push( vertStruct_1 );
+
+ //counts how often a node is saved in the list of lists
+ if(count_in_list == 1)
+ {
+ smiError =
+ SMI_Get_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndices[vertStruct_1.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n",
+ smiError);
+
+ value_q5 = value_q5 + 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, 1,
+ const_cast<DegreeOfFreedom*>(& nodeIndices[vertStruct_1.VertNum]),
+ &value_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n",
+ smiError);
+ }
+
+ if(value_q6-1 < Index)
+ {
+ Continue = true;
+ }
+
+ valuesINT[VertNum_2_in] = 1;
+
+ smiError =
+ SMI_Set_quantity_values(1, 1, 1, SMI_TYPE_INT, dim+1, 1,
+ &ElNum_in,
+ //&valuesINT);
+ valuesINT);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Set_quantity_values() failed with error %d\n", smiError);
+ }
+ }
+
+ delete [] valuesINT;
+ delete [] locVert;
+
+ return Continue;
+}
+
+bool
+HL_SignedDistLevels::getNextNeighbour (const int ElNum_in,
+ const int Vert_Up_in,
+ const int VertNum_in,
+ int &ElNum_out,
+ int &VertNum_1_out,
+ int &VertNum_2_out)
+{
+ int *neighbour = new int[dim+1];
+ int *oppVertices = new int[dim+1];
+ int *nodeIndicesOfElem;
+ int smiError;
+
+ smiAdapter->getNeighbourInfo(ElNum_in, neighbour,oppVertices);
+
+ ElNum_out = neighbour [VertNum_in];
+
+ if(ElNum_out != -1)
+ {
+ smiError =
+ SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&(ElNum_out)),
+ NULL, &nodeIndicesOfElem, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ VertNum_1_out = oppVertices[VertNum_in];
+ for (int i=0; i<=dim; i++)
+ {
+ if (nodeIndicesOfElem[i] != Vert_Up_in && i != VertNum_1_out)
+ {
+ VertNum_2_out = i;
+ }
+ }
+ }
+
+ delete [] neighbour;
+ delete [] oppVertices;
+
+ if (ElNum_out == -1)
+ {
+ return true;
+ }
+
+ return false;
+
+}
+
+int
+HL_SignedDistLevels::getNext_node_l_r (int elem_l_r_in,
+ int neighbour_l_r_in,
+ int node_l_r_in,
+ const int Vert)
+{
+ FUNCNAME("getNext_node_l_r()");
+
+ int *nodeIndicesOfElem;
+ int *nodeIndicesOfNeighbour = new int[dim+1];
+ int globalNextNode_l_r_out;
+
+ int smiError =
+ SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&neighbour_l_r_in),
+ NULL, &nodeIndicesOfElem, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ for(int i=0; i<=dim; i++)
+ {
+ nodeIndicesOfNeighbour[i] = nodeIndicesOfElem[i];
+ }
+
+ smiError =
+ SMI_Get_elems(1, 1, 1, const_cast<DegreeOfFreedom*>(&elem_l_r_in),
+ NULL, &nodeIndicesOfElem, NULL, NULL);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_elems() failed with error %d\n", smiError);
+
+ for (int i=0; i<=dim; i++)
+ {
+ if(nodeIndicesOfElem[i] != Vert && i != node_l_r_in)
+ {
+ globalNextNode_l_r_out = nodeIndicesOfElem[i];
+ }
+ }
+ for (int i=0; i<=dim; i++)
+ {
+ if(nodeIndicesOfNeighbour[i] == globalNextNode_l_r_out)
+ {
+ return i;
+ }
+ }
+
+ ERROR_EXIT("should never be reached !\n");
+ return 0;
+}
+
+void
+HL_SignedDistLevels::print_quantity_5 (int cntr)
+{
+ int numNodes;
+ int*nodeIndices;
+ double *coords = new double[dim];
+
+ int smiError =
+ SMI_Get_all_nodes(1, 1, const_cast<DegreeOfFreedom*>(&numNodes),
+ &nodeIndices);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_alls() failed with error %d\n", smiError);
+
+ int *values_q5 = new int [numNodes];
+
+ std::string q5_OutFile_0;
+ std::string q5_OutFile_1;
+ std::string q5_OutFile_2;
+ std::string q5_OutFile_3;
+ std::string q5_OutFile_4;
+
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_0", &q5_OutFile_0);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_1", &q5_OutFile_1);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_2", &q5_OutFile_2);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_3", &q5_OutFile_3);
+ GET_PARAMETER(0, "SignedDist->count saving->output->filename_4", &q5_OutFile_4);
+
+
+ char cntrStr[20];
+ sprintf(cntrStr, "%d", cntr);
+ q5_OutFile_0 += cntrStr;
+ q5_OutFile_1 += cntrStr;
+ q5_OutFile_2 += cntrStr;
+ q5_OutFile_3 += cntrStr;
+ q5_OutFile_4 += cntrStr;
+
+ cout << "count saving Ausgabe-Datei_0: " << q5_OutFile_0.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_1: " << q5_OutFile_1.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_2: " << q5_OutFile_2.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_3: " << q5_OutFile_3.c_str() << "\n\n";
+ cout << "count saving Ausgabe-Datei_4: " << q5_OutFile_4.c_str() << "\n\n";
+
+ ofstream out_0 (q5_OutFile_0.c_str());
+ ofstream out_1 (q5_OutFile_1.c_str());
+ ofstream out_2 (q5_OutFile_2.c_str());
+ ofstream out_3 (q5_OutFile_3.c_str());
+ ofstream out_4 (q5_OutFile_4.c_str());
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 5, SMI_TYPE_INT, 1, numNodes, nodeIndices,
+ values_q5);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ for (int i=0; i<numNodes; i++)
+ {
+ smiError =
+ SMI_Get_nodes (1, 1, 1, dim,
+ const_cast<DegreeOfFreedom*>(&nodeIndices[i]),
+ coords);
+
+ if(values_q5[i] == 0)
+ {
+ out_0<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] == 1)
+ {
+ out_1<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] == 2)
+ {
+ out_2<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] == 3)
+ {
+ out_3<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+ if(values_q5[i] >= 4)
+ {
+ out_4<<coords[0]<<" "<<coords[1]<<"\n";
+ }
+
+
+ }
+
+ out_0.close();
+ out_1.close();
+ out_2.close();
+ out_3.close();
+ out_4.close();
+
+ delete [] coords;
+ delete [] values_q5;
+}
+
+void
+HL_SignedDistLevels::print_quantity_6 ()
+{
+ std::string smiOutFile2;
+ GET_PARAMETER(0, "SignedDist->print levels->filename2", &smiOutFile2);
+ cout << "count levels Ausgabe-Datei: " << smiOutFile2.c_str() << "\n\n";
+
+ ofstream out2 (smiOutFile2.c_str());
+
+ int numElems;
+ int *elemIndices;
+
+ int smiError = SMI_Get_all_elems(1, 1, &numElems, &elemIndices);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_all_elems() failed with error %d\n", smiError);
+
+ int *value_q6 = new int [numElems];
+
+ smiError =
+ SMI_Get_quantity_values(1, 1, 6, SMI_TYPE_INT, 1, numElems, elemIndices,
+ value_q6);
+ TEST_EXIT(smiError == SMI_OK)
+ ("SMI_Get_quantity_values() failed with error %d\n", smiError);
+
+ int cl=0;
+ for (int i=0; i<numElems; i++)
+ {
+ out2<<elemIndices[i]<<" "<<value_q6[i]<<"\n";
+
+ if(value_q6[i] == 100000)
+ {
+ cl=cl +1;
+ }
+ }
+
+ out2.close();
+}
diff --git a/AMDiS/Reinit/src/HL_SignedDistLevels.h b/AMDiS/Reinit/src/HL_SignedDistLevels.h
new file mode 100644
index 00000000..aa02a31d
--- /dev/null
+++ b/AMDiS/Reinit/src/HL_SignedDistLevels.h
@@ -0,0 +1,234 @@
+#ifndef HL_SIGNEDDISTLEVELS
+#define HL_SIGNEDDISTLEVELS
+
+#include "ElInfo.h"
+#include "FixVec.h"
+// #include "MemoryManager.h"
+#include "Traverse.h"
+#include <queue>
+#include <vector>
+#include <time.h>
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+#include "ElementUpdate.h"
+#include "ElementUpdate_2d.h"
+#include "ElementUpdate_3d.h"
+#include "HL_SignedDist.h"
+#include "VelocityExt.h"
+
+#include "SMIAdapter.h"
+#include "smi.h"
+
+using namespace AMDiS;
+
+typedef struct
+{
+ int ElNum;
+ int VertNum;
+} Vert_Struct;
+
+class HL_SignedDistLevels : public HL_SignedDist
+{
+ public:
+ MEMORY_MANAGED(HL_SignedDistLevels);
+
+ /**
+ * Constructor.
+ */
+ HL_SignedDistLevels(const char *name_,
+ int dim_,
+ bool doVelocityExt = false,
+ Flag velExtType_ = VEL_EXT)
+ : HL_SignedDist(name_, dim_, doVelocityExt, velExtType_),
+ smiAdapter(NULL)
+ {
+ FUNCNAME("HL_SignedDistLevels::HL_SignedDistLevels");
+
+ // ===== Read parameters from init file. =====
+ GET_PARAMETER(0,name + "->tolerance", "%f", &tol);
+ GET_PARAMETER(0,name + "->count_how_often_saved_in_list", "%d", &count_in_list);
+ GET_PARAMETER(0,name + "->save_in_list->the ..th", "%d", &print_in_list);
+ GET_PARAMETER(0,name + "->save_in_list->after the ..th traversing of the list", "%d", &print_in_list_2);
+ GET_PARAMETER(0,name + "->print_level", "%d", &chosen_level);
+ GET_PARAMETER(0,name + "->print_level_yes_no", "%d", &print_level);
+
+ TEST_EXIT(tol > 0)("illegal tolerance !\n");
+ };
+
+ protected:
+ /**
+ * Initializes the boundary: calculation of the distance of boundary
+ * vertices to the interface.
+ * Interface is given by \ref lS_DOF and result is stored in
+ * \ref sD_DOF.
+ */
+ void initializeBoundary();
+
+ /**
+ * function for transvering Mesh to SMi an Adding the quantities
+ **/
+ void Mesh_to_SMI_and_quantity ();
+
+ /**
+ * function for creating the first list, in which elementnumber of
+ * boundary-elements are saved
+ */
+ void createFirstList(const int elStatus,
+ const int *elVertStatusVec,
+ int ElNum,
+ const int NumVertIntPoints);
+
+ /**
+ * Calculates the distance function and stores result in sD_DOF.
+ * Requirement: The boundary values are already set in \ref sD_DOF.
+ */
+ void HL_updateIteration();
+
+ /**
+ * function for traversing the list "List_Element"
+ * returns 1 if min one element is set in the new list
+ * returns 0 if no element is set in the new list
+ */
+ int traverseListElement();
+
+ /**
+ * function for creating all levels higher than 1 an saving pairs of
+ * elements and node in the right list
+ **/
+ void createLevels ();
+
+ /**
+ * function collects the next neighbours, creates an according list
+ * and saves the pair in it
+ * returns true if one or more neighbours exist and false if not
+ **/
+ bool collectNeighbours_setLevels (const int Element,
+ const int currentIndex,
+ bool *elementInNewListSet);
+
+ /**
+ * function for traversing the list of list, named "Level_"
+ **/
+ void traversingListLevel_ ( DOFVector<double> *boundVal_DOF );
+
+ /**
+ * calls the next neighbours and puts them into the right list
+ * returns true if an element had to be saved in an lower level
+ * than the current index
+ * this means the loop over all lists Level_[i] has to be repaeated again
+ * "Vert" is a local node
+ **/
+ bool search_and_include_comb(int ElNum,
+ int Vert,
+ int *nodeIndices,
+ const int Index);
+
+ /**
+ * checking whether the element "ElNum_in" has to be included into
+ * the second list
+ * if yes it will be included
+ * returns true if the level of the element is smaler than he current index
+ * returns false if the level is greater thean the current index
+ * "Vert_1_in" and "Vert_2_in" are localvertizes
+ */
+ bool includeIntoList (int ElNum_in,
+ int VertNum_1_in,
+ int VertNum_2_in,
+ const int Index);
+
+ /**
+ * gets the neighbour according to the node "VertNum_in" and two of
+ * its nodes (in local koordinates)
+ * returns "false" if the called neighbour exists and "true" if not
+ * Vert_Up_in is given in global coordinates, the other points are
+ * given in local coordinates
+ */
+ bool getNextNeighbour (const int ElNum_in,
+ const int Vert_Up_in,
+ const int VertNum_in,
+ int &ElNum_out,
+ int &VertNum_1_out,
+ int &VertNum_2_out);
+
+ /**
+ * function needed in the function "traversingListELVert"
+ **/
+ int getNext_node_l_r (int elem_l_r_in,
+ int neighbour_l_r_in,
+ int node_l_r_in,
+ const int Vert);
+
+ //==================================
+ void print_quantity_5 (int cntr);
+
+ /**
+ * function for printinq the quantity 6
+ * attention: we need at least read_transaction in smi for calling
+ * this function
+ **/
+ void print_quantity_6 ();
+ //====================================
+
+ protected:
+ /**
+ * Tolerance for Hopf-Lax update iteration loop.
+ */
+ double tol;
+
+ /**
+ * is needed for transfering the mesh to SMI
+ */
+ SMIAdapter *smiAdapter;
+
+ /**
+ * in this list boundary-elements are saved;
+ * is needed for creating the list "List_El_Vert"
+ */
+ queue<int> List_Element;
+
+ /**
+ * in this list structs filled with element-number and vertex-number
+ * are saved;
+ * is needed for traversing the mesh efficiently
+ */
+ queue<Vert_Struct> Level;
+
+ queue<Vert_Struct> helpLevel;
+
+ queue<Vert_Struct> List_El_Vert; //Listen noch einmal kontrollieren!
+
+ /**
+ * in this vector the lists of different levels are saved
+ **/
+ vector<queue<Vert_Struct> > Level_;
+
+ /**
+ * 0 ->do not count updates
+ * 1 ->count updates
+ **/
+ int count_updates;
+
+ /**
+ * 1 -> it will be count how often a node is saved in the second list
+ * 0 ->it will not be count
+ **/
+ int count_in_list;
+
+ int print_in_list;
+ int print_in_list_2;
+
+ /**
+ * level which will be printed in a file
+ **/
+ int chosen_level;
+
+ /**
+ * 1->a special level will be printed in a file
+ * 0->no level is printed in a file
+ **/
+ int print_level;
+};
+
+#endif // HL_SIGNEDDISTLEVELS
diff --git a/AMDiS/Reinit/src/HL_SignedDistTraverse.cc b/AMDiS/Reinit/src/HL_SignedDistTraverse.cc
new file mode 100644
index 00000000..534db935
--- /dev/null
+++ b/AMDiS/Reinit/src/HL_SignedDistTraverse.cc
@@ -0,0 +1,275 @@
+#include "HL_SignedDistTraverse.h"
+
+#include "VelocityExtFromVelocityField.h"
+
+void
+HL_SignedDistTraverse::initializeBoundary()
+{
+ FUNCNAME("HL_SignedDistTraverse::initializeBoundary()");
+
+ // ===== All non-boundary vertices are initialized with "infinity". =====
+ sD_DOF->set(inftyValue);
+
+ // ===== Traverse mesh and initialize boundary elements. =====
+ TraverseStack stack;
+ FixVec<double, VERTEX> distVec(dim, NO_INIT);
+ int elStatus;
+
+ const int nBasFcts = feSpace->getBasisFcts()->getNumber();
+ DegreeOfFreedom *locInd = GET_MEMORY(DegreeOfFreedom, nBasFcts);
+
+ ElInfo *elInfo;
+ if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
+ elInfo = stack.traverseFirst(feSpace->getMesh(),
+ -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_GRD_LAMBDA);
+ }
+ else {
+ elInfo = stack.traverseFirst(feSpace->getMesh(),
+ -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_COORDS);
+ }
+ while(elInfo) {
+
+ // Set elInfo in case velocity extension from velocity field is used.
+ if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
+ ((VelocityExtFromVelocityField *)velExt)->setElInfo(elInfo);
+ }
+
+ // Get local indices of vertices.
+ feSpace->getBasisFcts()->getLocalIndices(
+ const_cast<Element *>(elInfo->getElement()),
+ const_cast<DOFAdmin *>(feSpace->getAdmin()),
+ locInd);
+
+ // Get element status.
+ elStatus = elLS->createElementLevelSet(elInfo);
+
+ // Is element cut by the interface ?
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+
+ // Reset element distance vector.
+ for (int i=0; i<=dim; ++i) {
+ distVec[i] = inftyValue;
+ }
+
+ // Mark all vertices as boundary vertices.
+ for (int i=0; i<=dim; ++i) {
+ (*bound_DOF)[locInd[i]] = 1.0;
+ }
+
+ // Calculate distance for all vertices.
+ if (bndElDist->calcDistOnBoundaryElement(elInfo, distVec) !=
+ ElementLevelSet::LEVEL_SET_BOUNDARY) {
+ ERROR_EXIT("error in distance calculation !\n");
+ }
+ else {
+
+ // If distance is smaller, correct to new distance.
+ for (int i=0; i<=dim; ++i) {
+
+ // --> for test purposes:
+ if (distVec[i] > 1000) {
+ cout << "\nElement: " << elInfo->getElement()->getIndex() << ", Knoten: " << i << " , keine Randwertinitialisierung !\n";
+ }
+ // --> end: for test purposes
+
+ if ((*sD_DOF)[locInd[i]] > distVec[i]) {
+ (*sD_DOF)[locInd[i]] = distVec[i];
+ //If Distance is corrected, calculate new velocity.
+ if(velExt != NULL)
+ {
+ velExt->calcVelocityBoundary(locInd, i);
+ }
+ }
+ }
+ }
+ } // end of: elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY
+
+ elInfo = stack.traverseNext(elInfo);
+ } // end of: mesh traverse
+
+ FREE_MEMORY(locInd, DegreeOfFreedom, nBasFcts);
+
+ return;
+}
+
+void
+HL_SignedDistTraverse::HL_updateIteration()
+{
+ // ===== Create DOF vector for the last iteration step. =====
+ if (sDOld_DOF)
+ DELETE sDOld_DOF;
+ sDOld_DOF = NEW DOFVector<double>(feSpace, "sDOld_DOF");
+ sDOld_DOF->copy(const_cast<DOFVector<double> &>(*sD_DOF));
+
+ // ===== Gauss-Seidel or Jacobi iteration ? =====
+ if (GaussSeidelFlag) {
+ update_DOF = sD_DOF;
+ }
+ else {
+ update_DOF = sDOld_DOF;
+ }
+
+ // ===== Iteration loop: proceed until tolerance is reached. =====
+ TraverseStack stack;
+ ElInfo *elInfo;
+ tol_reached = false;
+ int itCntr = 0;
+ while (!tol_reached && itCntr != maxIt) {
+
+ ++itCntr;
+ tol_reached = true;
+
+ // ===== Traverse mesh: perform Hopf-Lax element update on
+ // each element. =====
+ elInfo = stack.traverseFirst(feSpace->getMesh(), -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_BOUND |
+ Mesh::FILL_COORDS);
+ while(elInfo) {
+ HL_elementUpdate(elInfo);
+
+ elInfo = stack.traverseNext(elInfo);
+ }
+
+ // ===== Is tolerance reached ? =====
+ tol_reached = checkTol();
+
+ sDOld_DOF->copy(const_cast<DOFVector<double> &>(*sD_DOF));
+ }
+
+ cout << "\n\nCalculation of signed distance function via mesh traverse iteration:\n";
+ if (GaussSeidelFlag) {
+ cout << "\tGauss-Seidel iteration\n";
+ }
+ else {
+ cout << "\tJacobi iteration\n";
+ }
+ cout << "\tnumber of iterations needed: " << itCntr << "\n\n";
+
+ return;
+}
+
+void
+HL_SignedDistTraverse::HL_elementUpdate(ElInfo *elInfo)
+{
+ FUNCNAME("HL_SignedDistTraverse::HL_elementUpdate()");
+
+ double update;
+
+ // ===== Get global indices of vertices of element. =====
+ const int nBasFcts = feSpace->getBasisFcts()->getNumber();
+ DegreeOfFreedom *locInd = GET_MEMORY(DegreeOfFreedom, nBasFcts);
+
+ feSpace->getBasisFcts()->getLocalIndices(
+ const_cast<Element *>(elInfo->getElement()),
+ const_cast<DOFAdmin *>(feSpace->getAdmin()),
+ locInd);
+
+ // ===== Hopf-Lax element update for each vertex of element. =====
+ for (int i=0; i<=dim; ++i) {
+
+ // ===== Calculate update for non-boundary vertex. =====
+ if ((*bound_DOF)[locInd[i]] < 1.e-15) {
+ //save permutation of vertexes for calculation of the velocity
+ if(velExt != NULL)
+ {
+ velExt->setPermutation(i, 1);
+ }
+ update = calcElementUpdate(elInfo, i, locInd);
+ // ---> for test purposes: count number of calculated updates
+ ++calcUpdate_Cntr;
+ // ---> end: for test purposes
+
+ // Calculates minimum of all element updates on elements
+ // containing vertex i.
+ if (update < (*sD_DOF)[locInd[i]]) {
+ (*sD_DOF)[locInd[i]] = update;
+ ///If Distance is corrected, calculate new velocity.
+ if(velExt != NULL)
+ {
+ velExt->calcVelocity(locInd, i);
+ }
+ // ---> for test purposes: count number of calculated updates
+ ++setUpdate_Cntr;
+ // ---> end: for test purposes
+ }
+ }
+ }
+
+ FREE_MEMORY(locInd, DegreeOfFreedom, nBasFcts);
+}
+
+double
+HL_SignedDistTraverse::calcElementUpdate(ElInfo *elInfo,
+ int nXh,
+ const DegreeOfFreedom *locInd)
+{
+ double update;
+ FixVec<WorldVector<double> *, VERTEX> elVert(dim, NO_INIT);
+ FixVec<double, VERTEX> uhVal(dim, NO_INIT);
+
+ // ===== Get local indices of element vertices (except xh). =====
+ int nYh, nZh, nWh;
+
+ nYh = (nXh + 1) % (dim+1);
+ nZh = (nXh + 2) % (dim+1);
+ if (dim == 3) nWh = (nXh + 3) % (dim+1);
+
+ // ===== Get world coordinates of vertices of element and their values
+ // of uh.
+ // The coordinates of the vertex the update is calculated for
+ // are stored at the end of the vector elVert. =====
+ switch (dim) {
+ case 2: elVert[0] = &(elInfo->getCoord(nYh));
+ elVert[1] = &(elInfo->getCoord(nZh));
+ elVert[2] = &(elInfo->getCoord(nXh));
+
+ uhVal[0] = (*update_DOF)[locInd[nYh]];
+ uhVal[1] = (*update_DOF)[locInd[nZh]];
+ uhVal[2] = (*update_DOF)[locInd[nXh]];
+
+ break;
+ case 3: elVert[0] = &(elInfo->getCoord(nYh));
+ elVert[1] = &(elInfo->getCoord(nZh));
+ elVert[2] = &(elInfo->getCoord(nWh));
+ elVert[3] = &(elInfo->getCoord(nXh));
+
+ uhVal[0] = (*update_DOF)[locInd[nYh]];
+ uhVal[1] = (*update_DOF)[locInd[nZh]];
+ uhVal[2] = (*update_DOF)[locInd[nWh]];
+ uhVal[3] = (*update_DOF)[locInd[nXh]];
+
+ break;
+ default: ERROR_EXIT("illegal dimension !\n");
+ break;
+ }
+
+ // ===== Calculate Hopf-Lax element update for vertex. =====
+ update = elUpdate->calcElementUpdate(elVert, uhVal);
+
+ return update;
+}
+
+bool
+HL_SignedDistTraverse::checkTol()
+{
+ DOFVector<double>::Iterator it_sD(sD_DOF, USED_DOFS);
+ DOFVector<double>::Iterator it_sDOld(sDOld_DOF, USED_DOFS);
+
+ for(it_sD.reset(), it_sDOld.reset();
+ !it_sD.end();
+ ++it_sD, ++it_sDOld) {
+ if ((*it_sDOld)-(*it_sD) > tol || (*it_sDOld)-(*it_sD) < 0) {
+ return false;
+ }
+ }
+
+ return true;
+}
diff --git a/AMDiS/Reinit/src/HL_SignedDistTraverse.h b/AMDiS/Reinit/src/HL_SignedDistTraverse.h
new file mode 100644
index 00000000..2f41d015
--- /dev/null
+++ b/AMDiS/Reinit/src/HL_SignedDistTraverse.h
@@ -0,0 +1,157 @@
+#ifndef HL_SIGNEDDISTTRAVERSE
+#define HL_SIGNEDDISTTRAVERSE
+
+#include "ElInfo.h"
+#include "FixVec.h"
+//#include "MemoryManager.h"
+#include "Traverse.h"
+
+#include "ElementLevelSet.h"
+
+#include "BoundaryElementDist.h"
+#include "ElementUpdate.h"
+#include "ElementUpdate_2d.h"
+#include "ElementUpdate_3d.h"
+#include "HL_SignedDist.h"
+#include "VelocityExt.h"
+
+using namespace AMDiS;
+
+class HL_SignedDistTraverse : public HL_SignedDist
+{
+ public:
+ MEMORY_MANAGED(HL_SignedDistTraverse);
+
+ /**
+ * Constructor.
+ */
+ HL_SignedDistTraverse(const char *name_,
+ int dim_,
+ bool doVelocityExt = false,
+ Flag velExtType_ = VEL_EXT)
+ : HL_SignedDist(name_, dim_, doVelocityExt, velExtType_),
+ sDOld_DOF(NULL),
+ update_DOF(NULL),
+ tol_reached(false)
+ {
+ FUNCNAME("HL_SignedDistTraverse::HL_SignedDistTraverse");
+
+ // ===== Read parameters from init file. =====
+ GET_PARAMETER(0, name + "->tolerance", "%f", &tol);
+ GET_PARAMETER(0, name + "->maximal number of iteration steps",
+ "%d", &maxIt);
+ GET_PARAMETER(0, name + "->Gauss-Seidel iteration", "%d",
+ &GaussSeidelFlag);
+
+ TEST_EXIT(tol > 0)("illegal tolerance !\n");
+
+ // ---> for test purposes: initialization of counter variables
+ calcUpdate_Cntr = 0;
+ setUpdate_Cntr = 0;
+ // ---> end: for test purposes
+ };
+
+ /**
+ * Destructor.
+ */
+ ~HL_SignedDistTraverse()
+ {
+ if (sDOld_DOF)
+ DELETE sDOld_DOF;
+
+ // ---> for test purposes: print result of update counting
+ printUpdateCntr();
+ // ---> end: for test purposes
+ };
+
+ protected:
+ /**
+ * Initializes the boundary: calculation of the distance of boundary
+ * vertices to the interface.
+ * Interface is given by lS_DOF and result is stored in
+ * sD_DOF.
+ */
+ void initializeBoundary();
+
+ /**
+ * Calculates the distance function and stores result in sD_DOF.
+ * Requirement: The boundary values are already set in sD_DOF.
+ */
+ void HL_updateIteration();
+
+ /**
+ * Hopf-Lax element update on element elInfo.
+ */
+ void HL_elementUpdate(ElInfo *elInfo);
+
+ /**
+ * Calculate the update for the vertex xh of element elInfo.
+ */
+ double calcElementUpdate(ElInfo *elInfo,
+ int nXh,
+ const DegreeOfFreedom *locInd);
+
+ /**
+ * In iteration loop: checks whether tolerance is reached.
+ * true - tolerance is reached
+ * false - tolerance is not reached
+ */
+ bool checkTol();
+
+ // ---> for test purposes: print result of update counting
+ void printUpdateCntr()
+ {
+ cout << "\n";
+ cout << "\tUpdate statistic: \n";
+ cout << "\t-----------------\n";
+ cout << "\tcalculated updates: " << calcUpdate_Cntr << "\n";
+ cout << "\tset updates: " << setUpdate_Cntr << "\n";
+ cout << "\n";
+ };
+ // ---> end: for test purposes
+
+ protected:
+ /**
+ * DOF vector for the last iteration step.
+ * Used during Jacobi iteration and to check whether tolerance is reached.
+ */
+ DOFVector<double> *sDOld_DOF;
+
+ /**
+ * Pointer to the DOF vector holding the values of the last iteration
+ * step.
+ * Gauss-Seidel iteration - update_DOF == sD_DOF
+ * Jacobi iteration - update_DOF == sDOld_DOF
+ */
+ DOFVector<double> *update_DOF;
+
+ /**
+ * Tolerance for Hopf-Lax update iteration loop.
+ */
+ double tol;
+
+ /**
+ * Flag showing whether tolerance tol is reached.
+ */
+ bool tol_reached;
+
+ /**
+ * Maximal number of mesh iterations for Hopf-Lax update.
+ */
+ int maxIt;
+
+ /**
+ * Indicates whether Gauss-Seidel or Jacobi iteration is used.
+ * 0 - Jacobi
+ * !=0 - Gauss-Seidel
+ */
+ int GaussSeidelFlag;
+
+ // ---> for test purposes: variables to count calculated and set updates
+ int calcUpdate_Cntr;
+ int setUpdate_Cntr;
+ // ---> end: for test purposes
+
+};
+
+#endif // HL_SIGNEDDISTTRAVERSE
diff --git a/AMDiS/Reinit/src/NormEps.cc b/AMDiS/Reinit/src/NormEps.cc
new file mode 100644
index 00000000..93adae67
--- /dev/null
+++ b/AMDiS/Reinit/src/NormEps.cc
@@ -0,0 +1,29 @@
+#include "NormEps.h"
+
+double NormEps::eps = 0.0;
+
+double
+NormEps::calcNormEps(const WorldVector<double> &x)
+{
+ FUNCNAME("NormEps::calcNormEps()");
+
+ static double result;
+
+ TEST_EXIT(eps > 1.e-15)("illegal eps for norm regularization !\n");
+ result = x*x + eps*eps;
+
+ return(sqrt(result));
+}
+
+double
+NormEps::calcNormEps2(const WorldVector<double> &x)
+{
+ FUNCNAME("NormEps::calcNormEps2()");
+
+ static double result;
+
+ TEST_EXIT(eps > 1.e-15)("illegal eps for norm regularization !\n");
+ result = x*x + eps*eps;
+
+ return(result);
+}
diff --git a/AMDiS/Reinit/src/NormEps.h b/AMDiS/Reinit/src/NormEps.h
new file mode 100644
index 00000000..1e83953b
--- /dev/null
+++ b/AMDiS/Reinit/src/NormEps.h
@@ -0,0 +1,46 @@
+#ifndef NORMEPS_H
+#define NORMEPS_H
+
+#include "FixVec.h"
+#include "Global.h"
+// #include "MemoryManager.h"
+#include "Parameters.h"
+
+using namespace AMDiS;
+
+class NormEps
+{
+ public:
+ MEMORY_MANAGED(NormEps);
+
+ /**
+ * Calculates regularized norm of WorldVector x.
+ */
+ static double calcNormEps(const WorldVector<double> &x);
+
+ /**
+ * Calculates square of regularized norm of WorldVector x.
+ */
+ static double calcNormEps2(const WorldVector<double> &x);
+
+ /**
+ * Set regularization epsilon eps.
+ */
+ static inline void setEps() {
+ FUNCNAME("NormEps::setEps()");
+
+ eps = 0.0;
+ GET_PARAMETER(0, "NormEps->epsilon for norm regularization", "%f", &eps);
+
+ TEST_EXIT(eps > 1.e-15)
+ ("illegal eps for norm regularization !\n");
+ };
+
+ protected:
+ /**
+ * Epsilon for regularization.
+ */
+ static double eps;
+};
+
+#endif // NORMEPS_H
diff --git a/AMDiS/Reinit/src/VelocityExt.cc b/AMDiS/Reinit/src/VelocityExt.cc
new file mode 100644
index 00000000..d7ca70b2
--- /dev/null
+++ b/AMDiS/Reinit/src/VelocityExt.cc
@@ -0,0 +1,260 @@
+#include "VelocityExt.h"
+
+void
+VelocityExt::calcVelocityBoundary(DegreeOfFreedom *lokInd, const int indexV)
+{
+ double tempV;
+
+ for (int nV=0; nV<nVelDOFs; ++nV) {
+
+ tempV = 0.0;
+ for(int i=0; i<=dim; i++)
+ {
+ tempV += (*(origVelDOF[nV]))[lokInd[i]]*lamVec[indexV][i];
+ }
+ (*(velDOF[nV]))[lokInd[indexV]]=tempV;
+ }
+}
+
+void
+VelocityExt::calcVelocity(DegreeOfFreedom *lokInd, const int indexV)
+{
+ double tempV;
+ int tempIndex;
+ if(indexFace>=0) //if there were several updates over element faces,
+ //take the coordinates that belong to the element face
+ //with the shortest distance to the interface
+ {
+ tempIndex=indexFace;
+ }
+ else //right coordinates are in lamVec[0]
+ {
+ tempIndex=0;
+ }
+ indexFace=-1;
+
+// // // test
+// int zero = permutation[0];
+// int one = permutation[1];
+// int two = permutation[2];
+// // int three = permutation[3];
+// // // test
+
+ for (int nV=0; nV<nVelDOFs; ++nV) {
+
+ tempV = 0.0;
+ for(int i=0; i<=dim; i++)
+ {
+ tempV += (*(velDOF[nV]))[lokInd[permutation[i]]]*lamVec[tempIndex][i];
+ }
+ (*(velDOF[nV]))[lokInd[indexV]]=tempV;
+ }
+}
+
+void
+VelocityExt::setBarycentricCoords_2D_boundary(const double &l_0, const double &l_1, const double &l_2, const int index)
+{
+ lamVec[index][0]=l_0;
+ lamVec[index][1]=l_1;
+ lamVec[index][2]=l_2;
+}
+
+void
+VelocityExt::setBarycentricCoords_2D(const double &l_0, const double &l_1, const double &l_2)
+{
+ setBarycentricCoords_2D_boundary(l_0, l_1, l_2, 0);
+}
+
+void
+VelocityExt::calcBarycentricCoords_2D(const double &c_delta, const double &c_alpha, const double &norm_zhminusyh, const double &norm_xhminusyh)
+{
+ double s_delta=sqrt(1-c_delta*c_delta);
+ double c_delta_alpha=c_delta*c_alpha+sqrt((1-c_alpha*c_alpha)*(1-c_delta*c_delta));
+// double c_delta_alpha=sqrt(c_delta*c_alpha+sqrt((1-c_alpha*c_alpha)*(1-c_delta*c_delta)));
+ double s_delta_alpha=sqrt(1-c_delta_alpha*c_delta_alpha);
+ double c=(s_delta_alpha*norm_xhminusyh)/s_delta;
+ double lambda=c/norm_zhminusyh;
+ setBarycentricCoords_2D(1-lambda, lambda, 0);
+}
+
+void
+VelocityExt::setBarycentricCoords_3D_boundary(const double &l_0, const double &l_1, const double &l_2, const double &l_3, const int index)
+{
+ lamVec[index][0]=l_0;
+ lamVec[index][1]=l_1;
+ lamVec[index][2]=l_2;
+ lamVec[index][3]=l_3;
+}
+
+void
+VelocityExt::calcBarycentricCoords_3D_boundary(const DimVec<double> sp1, const DimVec<double> sp2, const double lambda, int i)
+{
+ DimVec<double> tempVec(dim, NO_INIT);
+ tempVec[0] = sp1[0]+lambda*(sp2[0]-sp1[0]);
+ tempVec[1] = sp1[1]+lambda*(sp2[1]-sp1[1]);
+ tempVec[2] = sp1[2]+lambda*(sp2[2]-sp1[2]);
+ tempVec[3] = sp1[3]+lambda*(sp2[3]-sp1[3]);
+ setBarycentricCoords_3D_boundary(tempVec[0], tempVec[1], tempVec[2], tempVec[3], i);
+}
+
+void
+VelocityExt::copyAndExpandFaceCoords_3D(int vertNum,int index)
+{
+ double t0, t1, t2;
+ t0=lamVec[0][0];
+ t1=lamVec[0][1];
+ t2=lamVec[0][2];
+ switch(vertNum)
+ {
+ case 0: //coordinates belong to element face 0
+ setBarycentricCoords_3D_boundary(0,t0,t1,t2,index);
+ break;
+
+ case 1: //coordinates belong to element face 1
+ setBarycentricCoords_3D_boundary(t0,0,t1,t2,index);
+ break;
+
+ case 2: //coordinates belong to element face 2
+ setBarycentricCoords_3D_boundary(t0,t1,0,t2,index);
+ break;
+
+ default:
+ break;
+ }
+
+}
+
+void
+VelocityExt::setBarycentricCoords_3D(const double &l_0, const double &l_1, const double &l_2, const double &l_3)
+{
+ indexFace = -1;
+ setBarycentricCoords_3D_boundary(l_0, l_1, l_2, l_3, 0);
+}
+
+void
+VelocityExt::setIndexFaceWithShortestDist(int indexFace_)
+{
+ indexFace=indexFace_;
+}
+
+void
+VelocityExt::setPermutation(int vertNum, int mTrav)
+{
+ switch(dim)
+ {
+ case 2:
+ switch(mTrav)
+ {
+ case 0: //level
+ setPermutation(vertNum);
+ if(vertNum==1)
+ {
+ swapVertices(0,1);
+ }
+ break;
+
+ case 1: //traverse
+ setPermutation(vertNum);
+ break;
+
+ default:
+ break;
+ }
+ break;
+
+ case 3:
+ switch(mTrav)
+ {
+ case 1: //traverse
+ setPermutation(vertNum);
+ break;
+
+ default:
+ break;
+ }
+ break;
+
+ default:
+ break;
+ }
+}
+
+void
+VelocityExt::setPermutation(int vertNum)
+{
+ switch(dim)
+ {
+ case 2:
+ switch(vertNum)
+ {
+ case 0: //update is for vertex 0
+ setPermutation_2D(1,2,0);
+ break;
+
+ case 1: //update is for vertex 1
+ setPermutation_2D(2,0,1);
+ break;
+
+ case 2: //update is for vertex 2
+ setPermutation_2D(0,1,2);
+ break;
+
+ default:
+ break;
+ }
+ break;
+
+ case 3:
+ switch(vertNum)
+ {
+ case 0: //update is for vertex 0
+ setPermutation_3D(1,2,3,0);
+ break;
+
+ case 1: //update is for vertex 1
+ setPermutation_3D(2,3,0,1);
+ break;
+
+ case 2: //update is for vertex 2
+ setPermutation_3D(3,0,1,2);
+ break;
+
+ case 3: //update is for vertex 3
+ setPermutation_3D(0,1,2,3);
+ break;
+
+ default:
+ break;
+ }
+ break;
+
+ default:
+ break;
+ }
+}
+
+void
+VelocityExt::setPermutation_2D(int i_0, int i_1, int i_2)
+{
+ permutation[0]=i_0;
+ permutation[1]=i_1;
+ permutation[2]=i_2;
+}
+
+
+void
+VelocityExt::setPermutation_3D(int i0, int i1, int i2, int i3)
+{
+ permutation[0]=i0;
+ permutation[1]=i1;
+ permutation[2]=i2;
+ permutation[3]=i3;
+}
+
+void
+VelocityExt::swapVertices(int i1, int i2)
+{
+ int temp=permutation[i1];
+ permutation[i1]=permutation[i2];
+ permutation[i2]=temp;
+}
diff --git a/AMDiS/Reinit/src/VelocityExt.h b/AMDiS/Reinit/src/VelocityExt.h
new file mode 100644
index 00000000..0260fca2
--- /dev/null
+++ b/AMDiS/Reinit/src/VelocityExt.h
@@ -0,0 +1,258 @@
+#ifndef VELOCITYEXT_H
+#define VELOCITYEXT_H
+
+#include "AdaptInfo.h"
+#include "DOFVector.h"
+#include "Element.h"
+#include "FileWriter.h"
+#include "FixVec.h"
+// #include "MemoryManager.h"
+
+using namespace AMDiS;
+
+class VelocityExt
+{
+ public:
+ MEMORY_MANAGED(VelocityExt);
+
+ /**
+ * Constructor.
+ */
+ VelocityExt(int dim_)
+ : dim(dim_),
+ nVelDOFs(0),
+ lamVec(dim_,dim_+1,NO_INIT),
+ permutation(dim_,NO_INIT)
+ {
+ indexFace = -1;
+ };
+
+ /**
+ * Destructor.
+ */
+ virtual ~VelocityExt()
+ {};
+
+ /**
+ * Print velDOF to file.
+ */
+ void printVelDOF(AdaptInfo *adaptInfo, int i = 0)
+ {
+ FUNCNAME("VelocityExt::printVelDOF()");
+
+ TEST_EXIT(i < velDOF.size())("illegal index !\n");
+
+ FileWriter *fileWriter = NEW FileWriter(
+ "VelocityExt->velocity output",
+ (velDOF[i])->getFESpace()->getMesh(),
+ const_cast<DOFVector<double> *>(velDOF[i]));
+ fileWriter->writeFiles(adaptInfo, false);
+
+ DELETE fileWriter;
+ };
+
+ /**
+ * Print origVelDOF to file.
+ */
+ void printOrigVelDOF(AdaptInfo *adaptInfo, int i = 0)
+ {
+ FUNCNAME("VelocityExt::printOrigVelDOF()");
+
+ TEST_EXIT(i < origVelDOF.size())("illegal index !\n");
+
+ FileWriter *fileWriter = NEW FileWriter(
+ "VelocityExt->interface velocity output",
+ (origVelDOF[i])->getFESpace()->getMesh(),
+ const_cast<DOFVector<double> *>(origVelDOF[i]));
+ fileWriter->writeFiles(adaptInfo, false);
+
+ DELETE fileWriter;
+ };
+
+ /**
+ * Set velocity (one velocity vector).
+ */
+ void setVelocity(DOFVector<double> *origVelDOF_,
+ DOFVector<double> *velDOF_)
+ {
+ FUNCNAME("VelocityExt::setVelocity()");
+
+ nVelDOFs = 1;
+
+ TEST_EXIT(origVelDOF_)("illegal velocity vector origVelDOF !\n");
+ TEST_EXIT(velDOF_)("illegal velocity vector velDOF !\n");
+ TEST_EXIT(origVelDOF_->getFESpace() == velDOF_->getFESpace())
+ ("different fe spaces !\n");
+
+ origVelDOF.clear();
+ velDOF.clear();
+
+ origVelDOF.push_back(origVelDOF_);
+ velDOF.push_back(velDOF_);
+ };
+
+ /**
+ * Set velocity (multiple velocity vectors).
+ */
+ void setVelocity(std::vector<DOFVector<double> *> &origVelDOF_,
+ std::vector<DOFVector<double> *> &velDOF_)
+ {
+ FUNCNAME("VelocityExt::setVelocity()");
+
+ nVelDOFs = origVelDOF_.size();
+
+ TEST_EXIT(nVelDOFs > 0)("illegal number of velocity vectors !\n");
+ TEST_EXIT(nVelDOFs == velDOF_.size())("different sizes !\n");
+
+ for (int i=0; i<nVelDOFs; ++i) {
+ TEST_EXIT(origVelDOF_[i])("illegal velocity vector origVelDOF !\n");
+ TEST_EXIT(velDOF_[i])("illegal velocity vector velDOF !\n");
+ TEST_EXIT((origVelDOF_[i])->getFESpace() == (velDOF_[i])->getFESpace())
+ ("different fe spaces !\n");
+ }
+
+ origVelDOF = origVelDOF_;
+ velDOF = velDOF_;
+ };
+
+ /**
+ * Calculates the velocity for a vertex on a boundary element.
+ * lokInd gives the lokal index of the element vertices.
+ * indexV is the index of the vertex(with respect to the numeration on the element),
+ * for that the velocity shall be calculated.
+ */
+ virtual void calcVelocityBoundary(DegreeOfFreedom *lokInd,
+ const int indexV);
+
+ /**
+ * Calculates the velocity for a vertex on a non-boundary element.
+ * lokInd gives the lokal index of the element vertices.
+ * indexV is the index of the vertex(with respect to the numeration on the element),
+ * for that the velocity shall be calculated.
+ */
+ void calcVelocity(DegreeOfFreedom *lokInd, const int indexV);
+
+ /**
+ * Sets barycentric coordinates for 2D for a vertex on a boundary element.
+ * l_0. l_1, l_2 are the coordinates.
+ * indexV is the index of the vertex(with respect to the numeration on the element),
+ * for that the coordinates shall be stored.
+ */
+ void setBarycentricCoords_2D_boundary(const double &l_0, const double &l_1, const double &l_2, const int indexV);
+
+ /**
+ * Sets barycentric coordinates for 2D for a vertex on a non-boundary element.
+ * l_0. l_1, l_2 are the coordinates.
+ */
+ void setBarycentricCoords_2D(const double &l_0, const double &l_1, const double &l_2);
+
+ /**
+ * Calculates the barycentric coordinates for 2D for a vertex on a non-boundary element.
+ */
+ void calcBarycentricCoords_2D(const double &c_delta, const double &c_alpha, const double &norm_zhminusyh, const double &norm_xhminusyh);
+
+ /**
+ * Sets barycentric coordinates for 3D for a vertex on a boundary element
+ * (in case the coordinates are already known).
+ * l_0, l_1, l-2, l_3 are the coordinates.
+ * indexV is the index of the vertex(with respect to the numeration on the element),
+ * for that the coordinates shall be stored.
+ */
+ void setBarycentricCoords_3D_boundary(const double &l_0, const double &l_1, const double &l_2, const double &l_3, const int indexV);
+
+ /**
+ * Calculates the barycentric coordinates for 3D for a boundary vertex.
+ */
+ void calcBarycentricCoords_3D_boundary(const DimVec<double> sp1, const DimVec<double> sp2, const double lambda, int i);
+
+ /**
+ * Expands 3 coordinates from a face update to 4 coordinates and stores them in lamVec[index].
+ * vertNum is the index of the element face the coordinates belong to(0,1,or 2).
+ */
+ void copyAndExpandFaceCoords_3D(int vertNum,int index);
+
+/**
+ * Sets barycentric coordinates for 3D for a vertex on a non-boundary element.
+ * l_0, l_1, l-2, l_3 are the coordinates.
+ */
+ void setBarycentricCoords_3D(const double &l_0, const double &l_1, const double &l_2, const double &l_3);
+
+ /**
+ * Stores the index of element face with the shortest distance to the interface so far.
+ * (with respect to the order of calculation)
+ */
+ void setIndexFaceWithShortestDist(int indexC);
+
+ /**
+ * Sets the permutation of the vertices.
+ * vertNum is the vertex number(with respect to the numeration on the element),
+ * for that the update is calculated.
+ * mTrav is an index for determining the way of mesh-traverse.
+ * mTrav=0:level
+ * mTrav=1:traverse
+ */
+ void setPermutation(int vertNum, int mTrav);
+
+ /**
+ * Sets the permutation.
+ */
+ void setPermutation(int vertNum);
+
+ /**
+ * Sets the permutation of the vertices in 2D.
+ */
+ void setPermutation_2D(int i_0, int i_1, int i_2);
+
+ /**
+ * Sets the permutation of the vertices in 3D.
+ */
+ void setPermutation_3D(int i_0, int i_1, int i_2, int i_3);
+
+ /**
+ * Swaps two vertices in the permutation.
+ */
+ void swapVertices(int i1, int i2);
+
+ protected:
+ /**
+ * Original velocity vector.
+ */
+ std::vector<DOFVector<double> *> origVelDOF;
+
+ /**
+ * Dimension of mesh.
+ */
+ int dim;
+
+ /**
+ * Dof vector with extended velocity.
+ */
+ std::vector<DOFVector<double> *> velDOF;
+
+ /**
+ * Number of velocity vectors to be extended.
+ */
+ int nVelDOFs;
+
+ /**
+ * Vector with barycentric coordinates of points from which the velocity is extended.
+ * Whole vector is needed for the calculation on the boundary elements.
+ * For non-boundary elements the first place is used.
+ * If the update is calculated with the element faces, the places two-dim are used.
+ */
+ VectorOfFixVecs<DimVec<double> > lamVec;
+
+ /**
+ * If the update is calculated with the element faces, indexFace stores the number of the element
+ * face with the shortest distance to the interface.
+ * Otherwise indexFace is -1;
+ */
+ int indexFace;
+
+ /**
+ * List with the permutation of the vertices.
+ */
+ DimVec<int> permutation;
+};
+
+#endif // VELOCITYEXT_H
diff --git a/AMDiS/Reinit/src/VelocityExtFromVelocityField.cc b/AMDiS/Reinit/src/VelocityExtFromVelocityField.cc
new file mode 100644
index 00000000..d90c25b4
--- /dev/null
+++ b/AMDiS/Reinit/src/VelocityExtFromVelocityField.cc
@@ -0,0 +1,42 @@
+#include "VelocityExtFromVelocityField.h"
+
+// double NormEps::eps = 0.0;
+
+void
+VelocityExtFromVelocityField::calcVelocityBoundary(DegreeOfFreedom *locInd,
+ const int indexV)
+{
+ FUNCNAME("VelocityExtFromVelocityField::calcVelocityBoundary()");
+
+ // ===== Calculate normal velocity in element vertices. =====
+
+ // Get gradient of lSFct on element.
+ DimVec<double> lambda(dim, DEFAULT_VALUE, 0.0);
+ lambda[0] = 1.0;
+ WorldVector<double> elGrd;
+
+ for (int i=0; i<=dim; ++i)
+ lSFctVal[i] = (*lSFct)[locInd[i]];
+
+ basFcts->evalGrdUh(lambda, elInfo->getGrdLambda(), lSFctVal, &elGrd);
+
+ // Calculate normal velocity.
+ double fac = 1.0 / NormEps::calcNormEps(elGrd);
+ double sP;
+ for (int i=0; i<=dim; ++i) {
+ sP = 0.0;
+ for (int j=0; j<dim; ++j) {
+ sP += (*velField[j])[locInd[i]] * elGrd[j];
+ }
+ elNormalVel[i] = fac * sP;
+ }
+
+
+ // ===== Extend normal velocity from interface to element vertices. =====
+ double tempV=0.0;
+ for(int i=0; i<=dim; i++)
+ {
+ tempV += elNormalVel[i]*lamVec[indexV][i];
+ }
+ (*(velDOF[0]))[locInd[indexV]]=tempV;
+}
diff --git a/AMDiS/Reinit/src/VelocityExtFromVelocityField.h b/AMDiS/Reinit/src/VelocityExtFromVelocityField.h
new file mode 100644
index 00000000..d606a982
--- /dev/null
+++ b/AMDiS/Reinit/src/VelocityExtFromVelocityField.h
@@ -0,0 +1,181 @@
+#ifndef VELOCITYEXTFROMVELOCITYFIELD_H
+#define VELOCITYEXTFROMVELOCITYFIELD_H
+
+#include "AdaptInfo.h"
+#include "BasisFunction.h"
+#include "DOFVector.h"
+#include "ElInfo.h"
+// #include "MemoryManager.h"
+
+#include "NormEps.h"
+#include "VelocityExt.h"
+
+using namespace AMDiS;
+
+
+/* #ifndef NORMEPS_H */
+/* #define NORMEPS_H */
+/* ///////////////////////////////////////////////////////////////////////////// */
+/* // c l a s s N o r m E p s // */
+/* ///////////////////////////////////////////////////////////////////////////// */
+/* class NormEps */
+/* { */
+/* public: */
+/* MEMORY_MANAGED(NormEps); */
+
+/* /\** */
+/* * Calculates regularized norm of WorldVector x. */
+/* *\/ */
+/* static double calcNormEps(const WorldVector<double> &x) */
+/* { */
+/* FUNCNAME("NormEps::calcNormEps()"); */
+
+/* static double result; */
+
+/* TEST_EXIT(eps > 1.e-15)("illegal eps for norm regularization !\n"); */
+/* result = x*x + eps*eps; */
+
+/* return(std::sqrt(result)); */
+/* }; */
+
+/* /\** */
+/* * Calculates square of regularized norm of WorldVector x. */
+/* *\/ */
+/* static double calcNormEps2(const WorldVector<double> &x) */
+/* { */
+/* FUNCNAME("NormEps::calcNormEps2()"); */
+
+/* static double result; */
+
+/* TEST_EXIT(eps > 1.e-15)("illegal eps for norm regularization !\n"); */
+/* result = x*x + eps*eps; */
+
+/* return(result); */
+/* }; */
+
+/* /\** */
+/* * Set regularization epsilon eps. */
+/* *\/ */
+/* static inline void setEps() { */
+/* FUNCNAME("NormEps::setEps()"); */
+
+/* eps = 0.0; */
+/* GET_PARAMETER(0, "NormEps->epsilon for norm regularization", "%f", &eps); */
+
+/* TEST_EXIT(eps > 1.e-15) */
+/* ("illegal eps for norm regularization !\n"); */
+/* }; */
+
+/* protected: */
+/* /\** */
+/* * Epsilon for regularization. */
+/* *\/ */
+/* static double eps; */
+/* }; */
+/* #endif // NORMEPS_H */
+
+/////////////////////////////////////////////////////////////////////////////
+// c l a s s V e l o c i t y E x t F r o m V e l o c i t y F i e l d //
+/////////////////////////////////////////////////////////////////////////////
+class VelocityExtFromVelocityField : public VelocityExt
+{
+ public:
+ MEMORY_MANAGED(VelocityExtFromVelocityField);
+
+ /**
+ * Constructor.
+ */
+ VelocityExtFromVelocityField(int dim_)
+ : VelocityExt(dim_),
+ lSFct(NULL),
+ elNormalVel(dim_, NO_INIT),
+ basFcts(NULL)
+ {
+ FUNCNAME("VelocityExtFromVelocityField::VelocityExtFromVelocityField()");
+
+ lSFctVal = GET_MEMORY(double, dim+1);
+
+ // ===== set epsilon for norm regularization =====
+ NormEps::setEps();
+ };
+
+ /**
+ * Destructor.
+ */
+ ~VelocityExtFromVelocityField()
+ {
+ FREE_MEMORY(lSFctVal, double, dim+1);
+ };
+
+ /**
+ * Set velocity field.
+ */
+ void setVelocityField(std::vector<DOFVector<double> *> &velField_,
+ const DOFVector<double> *lSFct_,
+ DOFVector<double> *velDOF_)
+ {
+ FUNCNAME("VelocityExtFromVelocityField::setVelocityField()");
+
+ nVelDOFs = 1;
+
+ velField = velField_;
+ lSFct = lSFct_;
+ velDOF.clear();
+ velDOF.push_back(velDOF_);
+ origVelDOF.clear();
+
+ TEST_EXIT(lSFct)("level set function not defined !\n");
+ TEST_EXIT((int)velField.size() == dim)("illegal velocity field !\n");
+ TEST_EXIT(lSFct->getFESpace() == velDOF_->getFESpace())
+ ("different feSpaces !\n");
+
+ basFcts = lSFct->getFESpace()->getBasisFcts();
+ };
+
+ /**
+ * Adaption of VelocityExt::calcVelocityBoundary().
+ * Used to initialize normal velocity at interface elements.
+ */
+ void calcVelocityBoundary(DegreeOfFreedom *locInd, const int indexV);
+
+ /**
+ * Sets elInfo.
+ */
+ inline void setElInfo(ElInfo *elInfo_) {
+ elInfo = elInfo_;
+ };
+
+ protected:
+ /**
+ * Velocity field at interface.
+ */
+ std::vector<DOFVector<double> *> velField;
+
+ /**
+ * Interface is zero level set of level set function lSFct.
+ */
+ const DOFVector<double> *lSFct;
+
+ /**
+ * Normal velocity on a single element. Used to set normal velocity
+ * in velDOF on interface elements.
+ */
+ DimVec<double> elNormalVel;
+
+ /**
+ * Values of level set function in vertices of element.
+ */
+ double *lSFctVal;
+
+ /**
+ * Basis functions.
+ */
+ const BasisFunction *basFcts;
+
+ /**
+ * ElInfo used in calcVelocityBoundary().
+ */
+ ElInfo *elInfo;
+};
+
+#endif // VELOCITYEXTFROMVELOCITYFIELD_H
--
GitLab