Added reinit library to AMDiS repository.

AMDiS/Reinit/Makefile

+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

AMDiS/Reinit/src/BoundaryElementDist.cc

+#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];
+  }
+}

AMDiS/Reinit/src/BoundaryElementDist.h

+#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

AMDiS/Reinit/src/BoundaryElementEdgeDist.cc

+#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;
+}

AMDiS/Reinit/src/BoundaryElementEdgeDist.h

+#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

AMDiS/Reinit/src/BoundaryElementLevelSetDist.cc

+#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;
+}

AMDiS/Reinit/src/BoundaryElementLevelSetDist.h

+#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

AMDiS/Reinit/src/BoundaryElementNormalDist.cc

+#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);
+}

AMDiS/Reinit/src/BoundaryElementNormalDist.h

+#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

AMDiS/Reinit/src/BoundaryElementTopDist.cc

+#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];
+		}