diff --git a/AMDiS/Reinit/src/ElementUpdate_2d.cc b/AMDiS/Reinit/src/ElementUpdate_2d.cc
index 92d76aeef09abe703f032cfbe658dc796da1b26d..4d5584f4cc9026b9af95999b525f330a5e962d7e 100644
--- a/AMDiS/Reinit/src/ElementUpdate_2d.cc
+++ b/AMDiS/Reinit/src/ElementUpdate_2d.cc
@@ -1,51 +1,47 @@
#include "ElementUpdate_2d.h"
-double
-ElementUpdate_2d::calcElementUpdate(
- const FixVec<WorldVector<double> *, VERTEX> &vert,
- FixVec<double, VERTEX> &uhVal)
+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]);
+ WorldVector<double> &v0 = *(vert[0]);
+ WorldVector<double> &v1 = *(vert[1]);
+ WorldVector<double> &v2 = *(vert[2]);
+
+ int dim = Global::getGeo(WORLD);
+ for (int i = 0; i < dim; i++) {
+ xhminusYh[i] = v2[i] - v0[i];
+ zhminusYh[i] = v1[i] - v0[i];
+ xhminusZh[i] = v2[i] - v1[i];
+ }
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 delta = (uhVal[1] - uhVal[0]) / norm_zhminusYh;
double sP = xhminusYh * zhminusYh;
- c_alpha = sP / (norm_xhminusYh * norm_zhminusYh);
+ double c_alpha = sP / (norm_xhminusYh * norm_zhminusYh);
sP = xhminusZh * zhminusYh;
- c_beta = -sP / (norm_xhminusZh * norm_zhminusYh);
+ double 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);
- }
+ 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);
- }
+ 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))) *
+ (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);
- }
+ 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
index 3902bb8add68332445218fc37516b28490d75366..4a85a3895e38f54b11b4234d25b5fa031870bace 100644
--- a/AMDiS/Reinit/src/ElementUpdate_2d.h
+++ b/AMDiS/Reinit/src/ElementUpdate_2d.h
@@ -14,12 +14,15 @@ public:
: ElementUpdate(velExt_)
{}
- /**
- * Realization of ElementUpdate::calcElementUpdate.
- * Calculates the Bornemann update for an element.
+ /** \brief
+ * Realization of ElementUpdate::calcElementUpdate. Calculates the Bornemann
+ * update for an element.
*/
double calcElementUpdate(const FixVec<WorldVector<double> *, VERTEX> &vert,
FixVec<double, VERTEX> &uhVal);
+
+private:
+ WorldVector<double> xhminusYh, zhminusYh, xhminusZh;
};
#endif // ELEMENTUPDATE_2D_H
diff --git a/AMDiS/Reinit/src/HL_SignedDistTraverse.cc b/AMDiS/Reinit/src/HL_SignedDistTraverse.cc
index ace9bb587c0810defbd8d3d775062ba17f6c1497..bc782dc01ae5659e01f4a90891a4e51ae8042cc3 100644
--- a/AMDiS/Reinit/src/HL_SignedDistTraverse.cc
+++ b/AMDiS/Reinit/src/HL_SignedDistTraverse.cc
@@ -14,7 +14,8 @@ void HL_SignedDistTraverse::initializeBoundary()
int elStatus;
const int nBasFcts = feSpace->getBasisFcts()->getNumber();
- DegreeOfFreedom *locInd = new DegreeOfFreedom[nBasFcts];
+ if (locInd.size() < nBasFcts)
+ locInd.resize(nBasFcts);
ElInfo *elInfo;
if (velExt && velExtType.isSet(VEL_EXT_FROM_VEL_FIELD)) {
@@ -31,18 +32,18 @@ void HL_SignedDistTraverse::initializeBoundary()
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);
- }
+ 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);
+ const_cast<Element*>(elInfo->getElement()),
+ const_cast<DOFAdmin*>(feSpace->getAdmin()),
+ &(locInd[0]));
// Get element status.
elStatus = elLS->createElementLevelSet(elInfo);
@@ -51,85 +52,73 @@ void HL_SignedDistTraverse::initializeBoundary()
if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
// Reset element distance vector.
- for (int i=0; i<=dim; ++i) {
- distVec[i] = inftyValue;
- }
+ 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;
- }
+ 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 {
+ } else {
// If distance is smaller, correct to new distance.
- for (int i=0; i<=dim; ++i) {
+ for (int i = 0; i <= dim; i++) {
// --> for test purposes:
- if (distVec[i] > 1000) {
- cout << "\nElement: " << elInfo->getElement()->getIndex() << ", Knoten: " << i << " , keine Randwertinitialisierung !\n";
- }
+ 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);
- }
+ if (velExt != NULL)
+ velExt->calcVelocityBoundary(&(locInd[0]), i);
}
}
}
} // end of: elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY
elInfo = stack.traverseNext(elInfo);
- } // end of: mesh traverse
-
- delete [] locInd;
+ } // end of: mesh traverse
}
-void
-HL_SignedDistTraverse::HL_updateIteration()
+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");
+ 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;
- }
+ 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) {
-
+ 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) {
+ // ===== 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);
}
@@ -140,44 +129,39 @@ HL_SignedDistTraverse::HL_updateIteration()
}
cout << "\n\nCalculation of signed distance function via mesh traverse iteration:\n";
- if (GaussSeidelFlag) {
- cout << "\tGauss-Seidel iteration\n";
- }
- else {
+
+ if (GaussSeidelFlag)
+ cout << "\tGauss-Seidel iteration\n";
+ else
cout << "\tJacobi iteration\n";
- }
- cout << "\tnumber of iterations needed: " << itCntr << "\n\n";
- return;
+ cout << "\tnumber of iterations needed: " << itCntr << "\n\n";
}
-void
-HL_SignedDistTraverse::HL_elementUpdate(ElInfo *elInfo)
+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 = new DegreeOfFreedom[nBasFcts];
+ if (locInd.size() < nBasFcts)
+ locInd.resize(nBasFcts);
feSpace->getBasisFcts()->getLocalIndices(
const_cast<Element *>(elInfo->getElement()),
const_cast<DOFAdmin *>(feSpace->getAdmin()),
- locInd);
+ &(locInd[0]));
// ===== Hopf-Lax element update for each vertex of element. =====
- for (int i=0; i<=dim; ++i) {
+ 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);
+ if (velExt != NULL)
+ velExt->setPermutation(i, 1);
+
+ double update = calcElementUpdate(elInfo, i, &(locInd[0]));
// ---> for test purposes: count number of calculated updates
++calcUpdate_Cntr;
// ---> end: for test purposes
@@ -186,42 +170,37 @@ HL_SignedDistTraverse::HL_elementUpdate(ElInfo *elInfo)
// 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);
- }
+
+ // If Distance is corrected, calculate new velocity.
+ if(velExt != NULL)
+ velExt->calcVelocity(&(locInd[0]), i);
+
// ---> for test purposes: count number of calculated updates
++setUpdate_Cntr;
// ---> end: for test purposes
}
}
- }
-
- delete [] locInd;
+ }
}
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);
+ int nWh = 0;
+ int nYh = (nXh + 1) % (dim+1);
+ int 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));
+ case 2:
+ elVert[0] = &(elInfo->getCoord(nYh));
elVert[1] = &(elInfo->getCoord(nZh));
elVert[2] = &(elInfo->getCoord(nXh));
@@ -230,7 +209,8 @@ double HL_SignedDistTraverse::calcElementUpdate(ElInfo *elInfo,
uhVal[2] = (*update_DOF)[locInd[nXh]];
break;
- case 3: elVert[0] = &(elInfo->getCoord(nYh));
+ case 3:
+ elVert[0] = &(elInfo->getCoord(nYh));
elVert[1] = &(elInfo->getCoord(nZh));
elVert[2] = &(elInfo->getCoord(nWh));
elVert[3] = &(elInfo->getCoord(nXh));
@@ -246,9 +226,7 @@ double HL_SignedDistTraverse::calcElementUpdate(ElInfo *elInfo,
}
// ===== Calculate Hopf-Lax element update for vertex. =====
- update = elUpdate->calcElementUpdate(elVert, uhVal);
-
- return update;
+ return elUpdate->calcElementUpdate(elVert, uhVal);
}
bool HL_SignedDistTraverse::checkTol()
@@ -256,13 +234,9 @@ 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) {
+ 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
index bb0fe7c443d7ca0d830fa349901cf0d735a9aabc..32e7b3f5f9717b1d68790d58d7f571d17985b04b 100644
--- a/AMDiS/Reinit/src/HL_SignedDistTraverse.h
+++ b/AMDiS/Reinit/src/HL_SignedDistTraverse.h
@@ -24,7 +24,9 @@ public:
: HL_SignedDist(name_, dim_, doVelocityExt, velExtType_),
sDOld_DOF(NULL),
update_DOF(NULL),
- tol_reached(false)
+ tol_reached(false),
+ elVert(dim_, NO_INIT),
+ uhVal(dim_, NO_INIT)
{
FUNCNAME("HL_SignedDistTraverse::HL_SignedDistTraverse");
@@ -46,7 +48,7 @@ public:
~HL_SignedDistTraverse()
{
if (sDOld_DOF)
- DELETE sDOld_DOF;
+ delete sDOld_DOF;
// ---> for test purposes: print result of update counting
printUpdateCntr();
@@ -68,14 +70,10 @@ public:
*/
void HL_updateIteration();
- /**
- * Hopf-Lax element update on element elInfo.
- */
+ /// Hopf-Lax element update on element elInfo.
void HL_elementUpdate(ElInfo *elInfo);
- /**
- * Calculate the update for the vertex xh of element elInfo.
- */
+ /// Calculate the update for the vertex xh of element elInfo.
double calcElementUpdate(ElInfo *elInfo,
int nXh,
const DegreeOfFreedom *locInd);
@@ -96,7 +94,7 @@ public:
cout << "\tcalculated updates: " << calcUpdate_Cntr << "\n";
cout << "\tset updates: " << setUpdate_Cntr << "\n";
cout << "\n";
- };
+ }
// ---> end: for test purposes
protected:
@@ -114,19 +112,13 @@ public:
*/
DOFVector<double> *update_DOF;
- /**
- * Tolerance for Hopf-Lax update iteration loop.
- */
+ /// Tolerance for Hopf-Lax update iteration loop.
double tol;
- /**
- * Flag showing whether tolerance tol is reached.
- */
+ /// Flag showing whether tolerance tol is reached.
bool tol_reached;
- /**
- * Maximal number of mesh iterations for Hopf-Lax update.
- */
+ /// Maximal number of mesh iterations for Hopf-Lax update.
int maxIt;
/**
@@ -136,6 +128,12 @@ public:
*/
int GaussSeidelFlag;
+ std::vector<DegreeOfFreedom> locInd;
+
+ FixVec<WorldVector<double> *, VERTEX> elVert;
+
+ FixVec<double, VERTEX> uhVal;
+
// ---> for test purposes: variables to count calculated and set updates
int calcUpdate_Cntr;
int setUpdate_Cntr;