Commit 18086950 authored by Thomas Witkowski's avatar Thomas Witkowski
Browse files

* compositeFEM source adapted to AMDiS source layout

* deleted all gcc warnings from compiling compositeFEM
parent d3b88dd1
#include "CFE_Integration.h" #include "CFE_Integration.h"
#include "Mesh.h" #include "Mesh.h"
#include "SurfaceQuadrature.h" #include "SurfaceQuadrature.h"
#include "Traverse.h" #include "Traverse.h"
#include "ScalableQuadrature.h" #include "ScalableQuadrature.h"
#include "SubElInfo.h" #include "SubElInfo.h"
#include "SubPolytope.h" #include "SubPolytope.h"
double namespace AMDiS {
CFE_Integration::integrate_onNegLs(ElementFunction<double> *f,
ElementLevelSet *elLS, double
int deg, CFE_Integration::integrate_onNegLs(ElementFunction<double> *f,
Quadrature *q) ElementLevelSet *elLS,
{ int deg,
FUNCNAME("CFE_Integration::integrate_onNegLs()"); Quadrature *q)
{
int dim = elLS->getDim(); FUNCNAME("CFE_Integration::integrate_onNegLs()");
double int_val = 0.0;
double el_int_val; int dim = elLS->getDim();
double subEl_int_val; double int_val = 0.0;
VectorOfFixVecs<DimVec<double> > *intersecPts; double el_int_val;
SubPolytope *subPolytope; double subEl_int_val;
int numIntersecPts; VectorOfFixVecs<DimVec<double> > *intersecPts;
int iq; SubPolytope *subPolytope;
double val; int numIntersecPts;
int numQuadPts; int iq;
int vertex_interior; double val;
ScalableQuadrature *loc_scalQuad; int numQuadPts;
int numScalQuadPts; int vertex_interior;
bool subPolIsExterior = false; ScalableQuadrature *loc_scalQuad;
int elStatus; int numScalQuadPts;
Mesh *mesh = elLS->getMesh(); bool subPolIsExterior = false;
int elStatus;
// ===== Get quadratures. ===== Mesh *mesh = elLS->getMesh();
if (!q) {
q = Quadrature::provideQuadrature(dim, deg); // ===== Get quadratures. =====
} if (!q) {
numQuadPts = q->getNumPoints(); q = Quadrature::provideQuadrature(dim, deg);
loc_scalQuad = NEW ScalableQuadrature(q); }
numScalQuadPts = loc_scalQuad->getNumPoints(); numQuadPts = q->getNumPoints();
loc_scalQuad = NEW ScalableQuadrature(q);
numScalQuadPts = loc_scalQuad->getNumPoints();
// ===== Traverse mesh and calculate integral on each element. ===== // ===== Traverse mesh and calculate integral on each element. =====
TraverseStack stack; TraverseStack stack;
ElInfo *loc_elInfo = stack.traverseFirst(mesh, ElInfo *loc_elInfo = stack.traverseFirst(mesh,
-1, -1,
Mesh::CALL_LEAF_EL | Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS | Mesh::FILL_COORDS |
Mesh::FILL_DET); Mesh::FILL_DET);
while(loc_elInfo) { while(loc_elInfo) {
el_int_val = 0.0; el_int_val = 0.0;
subEl_int_val = 0.0; subEl_int_val = 0.0;
subPolIsExterior = false; subPolIsExterior = false;
// Check whether current element is cut by the zero level set. // Check whether current element is cut by the zero level set.
elStatus = elLS->createElementLevelSet(loc_elInfo); elStatus = elLS->createElementLevelSet(loc_elInfo);
if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) { if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
// ------------------------------------------------------------------- // -------------------------------------------------------------------
// Element is cut by the zero level set. // Element is cut by the zero level set.
// ------------------------------------------------------------------- // -------------------------------------------------------------------
// Create subelements. // Create subelements.
intersecPts = elLS->getElIntersecPoints(); intersecPts = elLS->getElIntersecPoints();
numIntersecPts = elLS->getNumElIntersecPoints(); numIntersecPts = elLS->getNumElIntersecPoints();
if (dim == 1 || (dim == 3 && numIntersecPts == 4)) { if (dim == 1 || (dim == 3 && numIntersecPts == 4)) {
// ----------------------------------------------------------------- // -----------------------------------------------------------------
// Subelement(s) are inside the domain with negative level set // Subelement(s) are inside the domain with negative level set
// function value. // function value.
// Get vertex with negative level set function value. // Get vertex with negative level set function value.
for (int i=0; i<=dim; ++i) { for (int i=0; i<=dim; ++i) {
if (elLS->getElVertLevelSetVec(i) < 0) { if (elLS->getElVertLevelSetVec(i) < 0) {
vertex_interior = i; vertex_interior = i;
break; break;
}
} }
}
subPolytope = NEW SubPolytope(loc_elInfo, subPolytope = NEW SubPolytope(loc_elInfo,
intersecPts, intersecPts,
numIntersecPts, numIntersecPts,
vertex_interior); vertex_interior);
} }
else { else {
// ----------------------------------------------------------------- // -----------------------------------------------------------------
// Subelement may be inside the domain with negative level set // Subelement may be inside the domain with negative level set
// function value as well as inside the domain with positive // function value as well as inside the domain with positive
// function value. // function value.
// //
// Whether a subelement is in the domain with negative or positive // Whether a subelement is in the domain with negative or positive
// level set function values is checked by the level set function // level set function values is checked by the level set function
// value of the first vertex of the subelement. (The subelements // value of the first vertex of the subelement. (The subelements
// are created in such a way that this vertex always is a vertex // are created in such a way that this vertex always is a vertex
// of the element and not an intersection point. Thus the level set // of the element and not an intersection point. Thus the level set
// function value of this vertex really is unequal to zero.) // function value of this vertex really is unequal to zero.)
subPolytope = NEW SubPolytope(loc_elInfo, subPolytope = NEW SubPolytope(loc_elInfo,
intersecPts, intersecPts,
numIntersecPts, numIntersecPts,
0); 0);
if(elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)) == if(elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)) ==
ElementLevelSet::LEVEL_SET_EXTERIOR) { ElementLevelSet::LEVEL_SET_EXTERIOR) {
subPolIsExterior = true; subPolIsExterior = true;
}
} }
}
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR); elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
// Calculate integral on subelement(s). // Calculate integral on subelement(s).
f->setElInfo(loc_elInfo); f->setElInfo(loc_elInfo);
for (std::vector<SubElInfo *>::iterator it = for (std::vector<SubElInfo *>::iterator it =
subPolytope->getSubElementsBegin(); subPolytope->getSubElementsBegin();
it != subPolytope->getSubElementsEnd(); it != subPolytope->getSubElementsEnd();
it++) { it++) {
loc_scalQuad->scaleQuadrature(**it);
loc_scalQuad->scaleQuadrature(**it); for (val = iq = 0; iq < numScalQuadPts; ++iq) {
val += loc_scalQuad->getWeight(iq)*(*f)(loc_scalQuad->getLambda(iq));
}
el_int_val += fabs((*it)->getDet()) * val;
}
for (val = iq = 0; iq < numScalQuadPts; ++iq) { // -------------------------------------------------------------------
val += loc_scalQuad->getWeight(iq)*(*f)(loc_scalQuad->getLambda(iq)); // In case the subelement is in the domain with positive level set
// function values:
// Calculate the integral on the element part with negative
// level set function values by substracting the integral on the
// subelement from the integral on the complete element.
if (subPolIsExterior) {
subEl_int_val = el_int_val;
el_int_val = 0.0;
f->setElInfo(loc_elInfo);
for (val = iq = 0; iq < numQuadPts; ++iq) {
val += q->getWeight(iq)*(*f)(q->getLambda(iq));
}
el_int_val = loc_elInfo->calcDet()*val - subEl_int_val;
} }
el_int_val += fabs((*it)->getDet()) * val;
// Free data.
DELETE subPolytope;
} }
else if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
// ------------------------------------------------------------------- // -------------------------------------------------------------------
// In case the subelement is in the domain with positive level set // Element is in the domain with negative level set function values.
// function values: // -------------------------------------------------------------------
// Calculate the integral on the element part with negative
// level set function values by substracting the integral on the elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
// subelement from the integral on the complete element.
if (subPolIsExterior) {
subEl_int_val = el_int_val;
el_int_val = 0.0;
f->setElInfo(loc_elInfo); f->setElInfo(loc_elInfo);
for (val = iq = 0; iq < numQuadPts; ++iq) { for (val = iq = 0; iq < numQuadPts; ++iq) {
val += q->getWeight(iq)*(*f)(q->getLambda(iq)); val += q->getWeight(iq)*(*f)(q->getLambda(iq));
} }
el_int_val = loc_elInfo->calcDet()*val - subEl_int_val; el_int_val = loc_elInfo->calcDet() * val;
}
// Free data.
DELETE subPolytope;
}
else if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
// -------------------------------------------------------------------
// Element is in the domain with negative level set function values.
// -------------------------------------------------------------------
elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
f->setElInfo(loc_elInfo);
for (val = iq = 0; iq < numQuadPts; ++iq) {
val += q->getWeight(iq)*(*f)(q->getLambda(iq));
} }
el_int_val = loc_elInfo->calcDet() * val;
}
int_val += el_int_val;
loc_elInfo = stack.traverseNext(loc_elInfo);
} // end of: mesh traverse int_val += el_int_val;
loc_elInfo = stack.traverseNext(loc_elInfo);
DELETE loc_scalQuad; } // end of: mesh traverse
return int_val; DELETE loc_scalQuad;
}
double
CFE_Integration::integrate_onZeroLs(ElementFunction<double> *f,
ElementLevelSet *elLS,
int deg,
Quadrature *q)
{
FUNCNAME("CFE_Integration::integrate_onZeroLs()");
int dim = elLS->getDim();
double int_val = 0.0;
VectorOfFixVecs<DimVec<double> > *intersecPts;
int numIntersecPts;
int iq;
double val;
SurfaceQuadrature *surfQuad;
int numQuadPts;
VectorOfFixVecs<DimVec<double> > tmpPts(dim, dim, NO_INIT);
DimVec<double> tmpPt(dim, DEFAULT_VALUE, 0.0);
int elStatus;
Mesh *mesh = elLS->getMesh();
// ===== Define default points for surface quadrature. =====
for (int i=0; i<dim; ++i) {
tmpPt[i] = 1.0;
tmpPts[i] = tmpPt;
tmpPt[i] = 0.0;
}
// ===== Get quadratures. ===== return int_val;
if (!q) {
q = Quadrature::provideQuadrature(dim-1, deg);
} }
surfQuad = NEW SurfaceQuadrature(q, tmpPts);
numQuadPts = surfQuad->getNumPoints();
double
CFE_Integration::integrate_onZeroLs(ElementFunction<double> *f,
ElementLevelSet *elLS,
int deg,
Quadrature *q)
{
FUNCNAME("CFE_Integration::integrate_onZeroLs()");
int dim = elLS->getDim();
double int_val = 0.0;
VectorOfFixVecs<DimVec<double> > *intersecPts;
int numIntersecPts;
int iq;
double val;
SurfaceQuadrature *surfQuad;
int numQuadPts;
VectorOfFixVecs<DimVec<double> > tmpPts(dim, dim, NO_INIT);
DimVec<double> tmpPt(dim, DEFAULT_VALUE, 0.0);
int elStatus;
Mesh *mesh = elLS->getMesh();
// ===== Define default points for surface quadrature. =====
for (int i=0; i<dim; ++i) {
tmpPt[i] = 1.0;
tmpPts[i] = tmpPt;
tmpPt[i] = 0.0;
}
// ===== Traverse mesh and calculate integral on each element cut by // ===== Get quadratures. =====
// the zero level set. ===== if (!q) {
TraverseStack stack; q = Quadrature::provideQuadrature(dim-1, deg);
}
ElInfo *loc_elInfo = stack.traverseFirst(mesh, surfQuad = NEW SurfaceQuadrature(q, tmpPts);
-1, numQuadPts = surfQuad->getNumPoints();
Mesh::CALL_LEAF_EL |
Mesh::FILL_COORDS |
Mesh::FILL_DET);
while(loc_elInfo) {
// Check whether current element is cut by the zero level set.
elStatus = elLS->createElementLevelSet(loc_elInfo);
if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
// Get intersection points.
intersecPts = elLS->getElIntersecPoints();
numIntersecPts = elLS->getNumElIntersecPoints();
// Calculate surface integral on intersection plane. // ===== Traverse mesh and calculate integral on each element cut by
f->setElInfo(loc_elInfo); // the zero level set. =====
TraverseStack stack;
// Note: The vector *intersecPts has always MAX_INTERSECTION_POINTS ElInfo *loc_elInfo = stack.traverseFirst(mesh,
// entries. -1,
for (int i=0; i<dim; ++i) Mesh::CALL_LEAF_EL |
tmpPts[i] = (*intersecPts)[i]; Mesh::FILL_COORDS |
Mesh::FILL_DET);
while(loc_elInfo) {
surfQuad->scaleSurfaceQuadrature(tmpPts); // Check whether current element is cut by the zero level set.
elStatus = elLS->createElementLevelSet(loc_elInfo);
for (val = iq = 0; iq < numQuadPts; ++iq) { if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
val += surfQuad->getWeight(iq)*(*f)(surfQuad->getLambda(iq));
} // Get intersection points.
int_val += calcSurfaceDet(loc_elInfo, tmpPts) * val; intersecPts = elLS->getElIntersecPoints();
numIntersecPts = elLS->getNumElIntersecPoints();
if (dim == 3 && numIntersecPts == 4) { // Calculate surface integral on intersection plane.
f->setElInfo(loc_elInfo);
// ----------------------------------------------------------------- // Note: The vector *intersecPts has always MAX_INTERSECTION_POINTS
// Intersection plane must be divided into two simplices. Calculate // entries.
// the surface integral for the second simplex. for (int i=0; i<dim; ++i)
// tmpPts[i] = (*intersecPts)[i];
// Note: The intersection points S0, S1, S2, S3 are supposed to be
// alligned in such a manner that a line through S1 and S2
// divides the intersection plane.
tmpPts[0] = (*intersecPts)[3];
surfQuad->scaleSurfaceQuadrature(tmpPts); surfQuad->scaleSurfaceQuadrature(tmpPts);
...@@ -262,33 +244,53 @@ CFE_Integration::integrate_onZeroLs(ElementFunction<double> *f, ...@@ -262,33 +244,53 @@ CFE_Integration::integrate_onZeroLs(ElementFunction<double> *f,
val += surfQuad->getWeight(iq)*(*f)(surfQuad->getLambda(iq)); val += surfQuad->getWeight(iq)*(*f)(surfQuad->getLambda(iq));
} }
int_val += calcSurfaceDet(loc_elInfo, tmpPts) * val; int_val += calcSurfaceDet(loc_elInfo, tmpPts) * val;
if (dim == 3 && numIntersecPts == 4) {
// -----------------------------------------------------------------
// Intersection plane must be divided into two simplices. Calculate
// the surface integral for the second simplex.
//
// Note: The intersection points S0, S1, S2, S3 are supposed to be
// alligned in such a manner that a line through S1 and S2
// divides the intersection plane.
tmpPts[0] = (*intersecPts)[3];
surfQuad->scaleSurfaceQuadrature(tmpPts);
for (val = iq = 0; iq < numQuadPts; ++iq) {
val += surfQuad->getWeight(iq)*(*f)(surfQuad->getLambda(iq));
}
int_val += calcSurfaceDet(loc_elInfo, tmpPts) * val;
}
} }
}
loc_elInfo = stack.traverseNext(loc_elInfo); loc_elInfo = stack.traverseNext(loc_elInfo);
} // end of: mesh traverse } // end of: mesh traverse
DELETE surfQuad; DELETE surfQuad;
return int_val; return int_val;
} }
double double
CFE_Integration::calcSurfaceDet(ElInfo *loc_elInfo, CFE_Integration::calcSurfaceDet(ElInfo *loc_elInfo,
VectorOfFixVecs<DimVec<double> > &surfVert) VectorOfFixVecs<DimVec<double> > &surfVert)
{ {
double surfDet; double surfDet;
int dim = surfVert[0].getSize()-1; int dim = surfVert[0].getSize()-1;
FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT); FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT);
// transform barycentric coords to world coords // transform barycentric coords to world coords
for(int i = 0; i < dim; i++) { for(int i = 0; i < dim; i++) {
loc_elInfo->coordToWorld(surfVert[i], &worldCoords[i]); loc_elInfo->coordToWorld(surfVert[i], &worldCoords[i]);
} }
// calculate determinant for surface // calculate determinant for surface
surfDet = ElInfo::calcDet(worldCoords); surfDet = ElInfo::calcDet(worldCoords);
return surfDet;
}
return surfDet;
} }
...@@ -4,40 +4,41 @@ ...@@ -4,40 +4,41 @@
#include "ElementFunction.h" #include "ElementFunction.h"
#include "MemoryManager.h" #include "MemoryManager.h"
#include "Quadrature.h" #include "Quadrature.h"
#include "ElementLevelSet.h" #include "ElementLevelSet.h"
using namespace AMDiS; namespace AMDiS {
class CFE_Integration
{
public:
MEMORY_MANAGED(CFE_Integration);
class CFE_Integration /**
{ * Calculates integral of function f on domain where level set function
public: * is negative.
MEMORY_MANAGED(CFE_Integration); */
static double integrate_onNegLs(ElementFunction<double> *f,
ElementLevelSet *elLS,
int deg = 1,
Quadrature *q = NULL);
/** /**
* Calculates integral of function f on domain where level set function * Calculates surface integral of function f on the zero level set.
* is negative. *
*/ * Note: Quadrature q is a quadrature formula for dimension dim-1.
static double integrate_onNegLs(ElementFunction<double> *f, */
ElementLevelSet *elLS, static double integrate_onZeroLs(ElementFunction<double> *f,
int deg = 1, ElementLevelSet *elLS,
Quadrature *q = NULL); int deg = 1,
Quadrature *q = NULL);
protected:
/**
* Calculates determinant for surface given through surfVert.
*/
static double calcSurfaceDet(ElInfo *loc_elInfo,
VectorOfFixVecs<DimVec<double> > &surfVert);
};
/**