Skip to content
Snippets Groups Projects
Commit 83c0388f authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

Some performance optimizations in reinit library.

parent 16cf3879
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
AMDiS/Reinit/src/ElementUpdate_2d.cc
+ 22
− 26
View file @ 83c0388f
#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;
......
AMDiS/Reinit/src/ElementUpdate_2d.h
+ 6
− 3
View file @ 83c0388f
......@@ -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
......
AMDiS/Reinit/src/HL_SignedDistTraverse.cc
+ 68
− 94
View file @ 83c0388f
......@@ -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;
}
AMDiS/Reinit/src/HL_SignedDistTraverse.h
+ 16
− 18
View file @ 83c0388f
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment