diff --git a/AMDiS/compositeFEM/src/CFE_Integration.cc b/AMDiS/compositeFEM/src/CFE_Integration.cc
index c320acf8dc1f91b580c817b364918559e909953c..8c4e0b29e6a2bfca8ecfef2ab0f406051c4d2083 100644
--- a/AMDiS/compositeFEM/src/CFE_Integration.cc
+++ b/AMDiS/compositeFEM/src/CFE_Integration.cc
@@ -1,260 +1,242 @@
#include "CFE_Integration.h"
-
#include "Mesh.h"
#include "SurfaceQuadrature.h"
#include "Traverse.h"
-
#include "ScalableQuadrature.h"
#include "SubElInfo.h"
#include "SubPolytope.h"
-double
-CFE_Integration::integrate_onNegLs(ElementFunction<double> *f,
- ElementLevelSet *elLS,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_Integration::integrate_onNegLs()");
-
- int dim = elLS->getDim();
- double int_val = 0.0;
- double el_int_val;
- double subEl_int_val;
- VectorOfFixVecs<DimVec<double> > *intersecPts;
- SubPolytope *subPolytope;
- int numIntersecPts;
- int iq;
- double val;
- int numQuadPts;
- int vertex_interior;
- ScalableQuadrature *loc_scalQuad;
- int numScalQuadPts;
- bool subPolIsExterior = false;
- int elStatus;
- Mesh *mesh = elLS->getMesh();
-
- // ===== Get quadratures. =====
- if (!q) {
- q = Quadrature::provideQuadrature(dim, deg);
- }
- numQuadPts = q->getNumPoints();
- loc_scalQuad = NEW ScalableQuadrature(q);
- numScalQuadPts = loc_scalQuad->getNumPoints();
+namespace AMDiS {
+
+ double
+ CFE_Integration::integrate_onNegLs(ElementFunction<double> *f,
+ ElementLevelSet *elLS,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_Integration::integrate_onNegLs()");
+
+ int dim = elLS->getDim();
+ double int_val = 0.0;
+ double el_int_val;
+ double subEl_int_val;
+ VectorOfFixVecs<DimVec<double> > *intersecPts;
+ SubPolytope *subPolytope;
+ int numIntersecPts;
+ int iq;
+ double val;
+ int numQuadPts;
+ int vertex_interior;
+ ScalableQuadrature *loc_scalQuad;
+ int numScalQuadPts;
+ bool subPolIsExterior = false;
+ int elStatus;
+ Mesh *mesh = elLS->getMesh();
+
+ // ===== Get quadratures. =====
+ if (!q) {
+ q = Quadrature::provideQuadrature(dim, deg);
+ }
+ numQuadPts = q->getNumPoints();
+ loc_scalQuad = NEW ScalableQuadrature(q);
+ numScalQuadPts = loc_scalQuad->getNumPoints();
- // ===== Traverse mesh and calculate integral on each element. =====
- TraverseStack stack;
+ // ===== Traverse mesh and calculate integral on each element. =====
+ TraverseStack stack;
- ElInfo *loc_elInfo = stack.traverseFirst(mesh,
- -1,
- Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET);
- while(loc_elInfo) {
- el_int_val = 0.0;
- subEl_int_val = 0.0;
- subPolIsExterior = false;
+ ElInfo *loc_elInfo = stack.traverseFirst(mesh,
+ -1,
+ Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET);
+ while(loc_elInfo) {
+ el_int_val = 0.0;
+ subEl_int_val = 0.0;
+ subPolIsExterior = false;
- // Check whether current element is cut by the zero level set.
- elStatus = elLS->createElementLevelSet(loc_elInfo);
+ // Check whether current element is cut by the zero level set.
+ 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.
- intersecPts = elLS->getElIntersecPoints();
- numIntersecPts = elLS->getNumElIntersecPoints();
+ // Create subelements.
+ intersecPts = elLS->getElIntersecPoints();
+ 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
- // function value.
+ // -----------------------------------------------------------------
+ // Subelement(s) are inside the domain with negative level set
+ // function value.
- // Get vertex with negative level set function value.
- for (int i=0; i<=dim; ++i) {
- if (elLS->getElVertLevelSetVec(i) < 0) {
- vertex_interior = i;
- break;
+ // Get vertex with negative level set function value.
+ for (int i=0; i<=dim; ++i) {
+ if (elLS->getElVertLevelSetVec(i) < 0) {
+ vertex_interior = i;
+ break;
+ }
}
- }
- subPolytope = NEW SubPolytope(loc_elInfo,
- intersecPts,
- numIntersecPts,
- vertex_interior);
- }
- else {
-
- // -----------------------------------------------------------------
- // Subelement may be inside the domain with negative level set
- // function value as well as inside the domain with positive
- // function value.
- //
- // Whether a subelement is in the domain with negative or positive
- // level set function values is checked by the level set function
- // value of the first vertex of the subelement. (The subelements
- // 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
- // function value of this vertex really is unequal to zero.)
-
- subPolytope = NEW SubPolytope(loc_elInfo,
- intersecPts,
- numIntersecPts,
- 0);
-
- if(elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)) ==
- ElementLevelSet::LEVEL_SET_EXTERIOR) {
- subPolIsExterior = true;
+ subPolytope = NEW SubPolytope(loc_elInfo,
+ intersecPts,
+ numIntersecPts,
+ vertex_interior);
+ }
+ else {
+
+ // -----------------------------------------------------------------
+ // Subelement may be inside the domain with negative level set
+ // function value as well as inside the domain with positive
+ // function value.
+ //
+ // Whether a subelement is in the domain with negative or positive
+ // level set function values is checked by the level set function
+ // value of the first vertex of the subelement. (The subelements
+ // 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
+ // function value of this vertex really is unequal to zero.)
+
+ subPolytope = NEW SubPolytope(loc_elInfo,
+ intersecPts,
+ numIntersecPts,
+ 0);
+
+ if(elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)) ==
+ ElementLevelSet::LEVEL_SET_EXTERIOR) {
+ subPolIsExterior = true;
+ }
}
- }
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
- // Calculate integral on subelement(s).
- f->setElInfo(loc_elInfo);
- for (std::vector<SubElInfo *>::iterator it =
- subPolytope->getSubElementsBegin();
- it != subPolytope->getSubElementsEnd();
- it++) {
+ // Calculate integral on subelement(s).
+ f->setElInfo(loc_elInfo);
+ for (std::vector<SubElInfo *>::iterator it =
+ subPolytope->getSubElementsBegin();
+ it != subPolytope->getSubElementsEnd();
+ 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
- // 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;
+ // -------------------------------------------------------------------
+ // 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 - subEl_int_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;
}
- 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;
-}
-
-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;
- }
+ DELETE loc_scalQuad;
- // ===== Get quadratures. =====
- if (!q) {
- q = Quadrature::provideQuadrature(dim-1, deg);
+ return int_val;
}
- 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
- // the zero level set. =====
- TraverseStack stack;
-
- ElInfo *loc_elInfo = stack.traverseFirst(mesh,
- -1,
- 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);
+ // ===== Get quadratures. =====
+ if (!q) {
+ q = Quadrature::provideQuadrature(dim-1, deg);
+ }
+ surfQuad = NEW SurfaceQuadrature(q, tmpPts);
+ numQuadPts = surfQuad->getNumPoints();
- if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
-
- // Get intersection points.
- intersecPts = elLS->getElIntersecPoints();
- numIntersecPts = elLS->getNumElIntersecPoints();
- // Calculate surface integral on intersection plane.
- f->setElInfo(loc_elInfo);
+ // ===== Traverse mesh and calculate integral on each element cut by
+ // the zero level set. =====
+ TraverseStack stack;
- // Note: The vector *intersecPts has always MAX_INTERSECTION_POINTS
- // entries.
- for (int i=0; i<dim; ++i)
- tmpPts[i] = (*intersecPts)[i];
+ ElInfo *loc_elInfo = stack.traverseFirst(mesh,
+ -1,
+ Mesh::CALL_LEAF_EL |
+ 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) {
- val += surfQuad->getWeight(iq)*(*f)(surfQuad->getLambda(iq));
- }
- int_val += calcSurfaceDet(loc_elInfo, tmpPts) * val;
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+
+ // Get intersection points.
+ intersecPts = elLS->getElIntersecPoints();
+ numIntersecPts = elLS->getNumElIntersecPoints();
- if (dim == 3 && numIntersecPts == 4) {
+ // Calculate surface integral on intersection plane.
+ f->setElInfo(loc_elInfo);
- // -----------------------------------------------------------------
- // 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];
+ // Note: The vector *intersecPts has always MAX_INTERSECTION_POINTS
+ // entries.
+ for (int i=0; i<dim; ++i)
+ tmpPts[i] = (*intersecPts)[i];
surfQuad->scaleSurfaceQuadrature(tmpPts);
@@ -262,33 +244,53 @@ CFE_Integration::integrate_onZeroLs(ElementFunction<double> *f,
val += surfQuad->getWeight(iq)*(*f)(surfQuad->getLambda(iq));
}
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
-CFE_Integration::calcSurfaceDet(ElInfo *loc_elInfo,
- VectorOfFixVecs<DimVec<double> > &surfVert)
-{
- double surfDet;
- int dim = surfVert[0].getSize()-1;
- FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT);
+ double
+ CFE_Integration::calcSurfaceDet(ElInfo *loc_elInfo,
+ VectorOfFixVecs<DimVec<double> > &surfVert)
+ {
+ double surfDet;
+ int dim = surfVert[0].getSize()-1;
+ FixVec<WorldVector<double>, VERTEX> worldCoords(dim-1, NO_INIT);
- // transform barycentric coords to world coords
- for(int i = 0; i < dim; i++) {
- loc_elInfo->coordToWorld(surfVert[i], &worldCoords[i]);
- }
+ // transform barycentric coords to world coords
+ for(int i = 0; i < dim; i++) {
+ loc_elInfo->coordToWorld(surfVert[i], &worldCoords[i]);
+ }
- // calculate determinant for surface
- surfDet = ElInfo::calcDet(worldCoords);
+ // calculate determinant for surface
+ surfDet = ElInfo::calcDet(worldCoords);
+
+ return surfDet;
+ }
- return surfDet;
}
diff --git a/AMDiS/compositeFEM/src/CFE_Integration.h b/AMDiS/compositeFEM/src/CFE_Integration.h
index 4c1fa133c7602a8f11c3f2d8fe35aa4e15d7abc2..3e2b22fa7ff8cdddf3c62b300af9eec89e52e37d 100644
--- a/AMDiS/compositeFEM/src/CFE_Integration.h
+++ b/AMDiS/compositeFEM/src/CFE_Integration.h
@@ -4,40 +4,41 @@
#include "ElementFunction.h"
#include "MemoryManager.h"
#include "Quadrature.h"
-
#include "ElementLevelSet.h"
-using namespace AMDiS;
+namespace AMDiS {
+
+ class CFE_Integration
+ {
+ public:
+ MEMORY_MANAGED(CFE_Integration);
-class CFE_Integration
-{
- public:
- MEMORY_MANAGED(CFE_Integration);
+ /**
+ * Calculates integral of function f on domain where level set function
+ * is negative.
+ */
+ 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
- * is negative.
- */
- static double integrate_onNegLs(ElementFunction<double> *f,
- ElementLevelSet *elLS,
- int deg = 1,
- Quadrature *q = NULL);
+ /**
+ * Calculates surface integral of function f on the zero level set.
+ *
+ * Note: Quadrature q is a quadrature formula for dimension dim-1.
+ */
+ static double integrate_onZeroLs(ElementFunction<double> *f,
+ ElementLevelSet *elLS,
+ int deg = 1,
+ Quadrature *q = NULL);
+ protected:
+ /**
+ * Calculates determinant for surface given through surfVert.
+ */
+ static double calcSurfaceDet(ElInfo *loc_elInfo,
+ VectorOfFixVecs<DimVec<double> > &surfVert);
+ };
- /**
- * Calculates surface integral of function f on the zero level set.
- *
- * Note: Quadrature q is a quadrature formula for dimension dim-1.
- */
- static double integrate_onZeroLs(ElementFunction<double> *f,
- ElementLevelSet *elLS,
- int deg = 1,
- Quadrature *q = NULL);
- protected:
- /**
- * Calculates determinant for surface given through surfVert.
- */
- static double calcSurfaceDet(ElInfo *loc_elInfo,
- VectorOfFixVecs<DimVec<double> > &surfVert);
-};
+}
#endif // AMDIS_CFE_INTEGRATION_H
diff --git a/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc b/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc
index 0cee8c39370f83ff495bbc90e79e0fd1ccb710b6..6ac83b9a1969055e9b74130609428d4a369dbb08 100644
--- a/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc
+++ b/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.cc
@@ -1,443 +1,566 @@
-#include "CFE_NormAndErrorFcts.h"
-
#include <vector>
+#include "CFE_NormAndErrorFcts.h"
#include "Mesh.h"
#include "Traverse.h"
-
#include "SubElInfo.h"
-double CFE_NormAndErrorFcts::L2_err_abs = 0.0;
-double CFE_NormAndErrorFcts::L2_u_norm = 0.0;
-double CFE_NormAndErrorFcts::H1_err_abs = 0.0;
-double CFE_NormAndErrorFcts::H1_u_norm = 0.0;
-
-double
-ElementL1Norm_Analyt::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- const WorldVector<double> *worldCoordsAtQP;
-
- for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
- val += q->getWeight(iq) * fabs((*f)(*worldCoordsAtQP));
- }
- double nrm = det * val;
+namespace AMDiS {
- return nrm;
-}
+ double CFE_NormAndErrorFcts::L2_err_abs = 0.0;
+ double CFE_NormAndErrorFcts::L2_u_norm = 0.0;
+ double CFE_NormAndErrorFcts::H1_err_abs = 0.0;
+ double CFE_NormAndErrorFcts::H1_u_norm = 0.0;
+
+ double
+ ElementL1Norm_Analyt::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ const WorldVector<double> *worldCoordsAtQP;
+
+ for (int iq = 0; iq < nQPts; ++iq) {
+ worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
+ val += q->getWeight(iq) * fabs((*f)(*worldCoordsAtQP));
+ }
+ double nrm = det * val;
-double
-ElementL2Norm_Analyt::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- const WorldVector<double> *worldCoordsAtQP;
-
- for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
- val += q->getWeight(iq) * sqr((*f)(*worldCoordsAtQP));
+ return nrm;
}
- double nrm = det * val;
- return nrm;
-}
+ double
+ ElementL2Norm_Analyt::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ const WorldVector<double> *worldCoordsAtQP;
+
+ for (int iq = 0; iq < nQPts; ++iq) {
+ worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
+ val += q->getWeight(iq) * sqr((*f)(*worldCoordsAtQP));
+ }
+ double nrm = det * val;
+
+ return nrm;
+ }
-double
-ElementH1Norm_Analyt::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- double norm_grd2;
- const WorldVector<double> *worldCoordsAtQP;
-
- for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
+ double
+ ElementH1Norm_Analyt::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ double norm_grd2;
+ const WorldVector<double> *worldCoordsAtQP;
+
+ for (int iq = 0; iq < nQPts; ++iq) {
+ worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
- norm_grd2 = 0.0;
- for (int j = 0; j < dim; j++)
- norm_grd2 += sqr(((*grd)(*worldCoordsAtQP))[j]);
+ norm_grd2 = 0.0;
+ for (int j = 0; j < dim; j++)
+ norm_grd2 += sqr(((*grd)(*worldCoordsAtQP))[j]);
- val += q->getWeight(iq) * norm_grd2;
+ val += q->getWeight(iq) * norm_grd2;
+ }
+ double nrm = det * val;
+
+ return nrm;
}
- double nrm = det * val;
- return nrm;
-}
+ double
+ ElementL1Norm_DOF::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ const double *dofAtQPs = dofVec->getVecAtQPs(elInfo,
+ q,
+ NULL,
+ NULL);
+
+ for (int iq = 0; iq < nQPts; ++iq) {
+ val += q->getWeight(iq) * fabs(dofAtQPs[iq]);
+ }
+ double nrm = det * val;
-double
-ElementL1Norm_DOF::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- const double *dofAtQPs = dofVec->getVecAtQPs(elInfo,
- q,
- NULL,
- NULL);
-
- for (int iq = 0; iq < nQPts; ++iq) {
- val += q->getWeight(iq) * fabs(dofAtQPs[iq]);
+ return nrm;
}
- double nrm = det * val;
- return nrm;
-}
+ double
+ ElementL2Norm_DOF::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ const double *dofAtQPs = dofVec->getVecAtQPs(elInfo,
+ q,
+ NULL,
+ NULL);
+
+ for (int iq = 0; iq < nQPts; ++iq) {
+ val += q->getWeight(iq) * sqr(dofAtQPs[iq]);
+ }
+ double nrm = det * val;
-double
-ElementL2Norm_DOF::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- const double *dofAtQPs = dofVec->getVecAtQPs(elInfo,
- q,
- NULL,
- NULL);
-
- for (int iq = 0; iq < nQPts; ++iq) {
- val += q->getWeight(iq) * sqr(dofAtQPs[iq]);
+ return nrm;
}
- double nrm = det * val;
- return nrm;
-}
-
-double
-ElementH1Norm_DOF::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- double norm_grd2;
- const WorldVector<double> *grdDofAtQPs = dofVec->getGrdAtQPs(elInfo,
- q,
- NULL,
- NULL);
-
- for (int iq = 0; iq < nQPts; ++iq) {
+ double
+ ElementH1Norm_DOF::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ double norm_grd2;
+ const WorldVector<double> *grdDofAtQPs = dofVec->getGrdAtQPs(elInfo,
+ q,
+ NULL,
+ NULL);
+
+ for (int iq = 0; iq < nQPts; ++iq) {
- norm_grd2 = 0.0;
- for (int j = 0; j < dim; ++j)
- norm_grd2 += sqr(grdDofAtQPs[iq][j]);
+ norm_grd2 = 0.0;
+ for (int j = 0; j < dim; ++j)
+ norm_grd2 += sqr(grdDofAtQPs[iq][j]);
- val += q->getWeight(iq) * norm_grd2;
- }
- double nrm = det * val;
+ val += q->getWeight(iq) * norm_grd2;
+ }
+ double nrm = det * val;
- return nrm;
-}
+ return nrm;
+ }
-double
-ElementL2Err::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- double val_nrm = 0.0;
- const double *uhAtQPs = uh->getVecAtQPs(elInfo,
- q,
- NULL,
- NULL);
- const WorldVector<double> *worldCoordsAtQP;
-
- for (int iq = 0; iq < nQPts; ++iq) {
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
- val += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP) - uhAtQPs[iq]);
+ double
+ ElementL2Err::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ double val_nrm = 0.0;
+ const double *uhAtQPs = uh->getVecAtQPs(elInfo,
+ q,
+ NULL,
+ NULL);
+ const WorldVector<double> *worldCoordsAtQP;
+
+ for (int iq = 0; iq < nQPts; ++iq) {
+ worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
+ val += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP) - uhAtQPs[iq]);
- if (relErr)
- val_nrm += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP));
- }
+ if (relErr)
+ val_nrm += q->getWeight(iq) * sqr((*u)(*worldCoordsAtQP));
+ }
- double err = det * val;
+ double err = det * val;
- if (relErr)
- nrmU += fac * det * val_nrm;
+ if (relErr)
+ nrmU += fac * det * val_nrm;
- return err;
-}
+ return err;
+ }
-double
-ElementH1Err::calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac)
-{
- double val = 0.0;
- double val_nrm = 0.0;
- double norm_err_grd2;
- double norm_grd2;
- const WorldVector<double> *grdUhAtQPs = uh->getGrdAtQPs(elInfo,
- q,
- NULL,
- NULL);
- const WorldVector<double> *worldCoordsAtQP;
-
- for (int iq = 0; iq < nQPts; ++iq) {
-
- worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
-
- norm_err_grd2 = 0.0;
- for (int j = 0; j < dim; ++j)
- norm_err_grd2 +=
- sqr(((*grdu)(*worldCoordsAtQP))[j] - grdUhAtQPs[iq][j]);
+ double
+ ElementH1Err::calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac)
+ {
+ double val = 0.0;
+ double val_nrm = 0.0;
+ double norm_err_grd2;
+ double norm_grd2;
+ const WorldVector<double> *grdUhAtQPs = uh->getGrdAtQPs(elInfo,
+ q,
+ NULL,
+ NULL);
+ const WorldVector<double> *worldCoordsAtQP;
+
+ for (int iq = 0; iq < nQPts; ++iq) {
+
+ worldCoordsAtQP = elInfo->coordToWorld(q->getLambda(iq), NULL);
+
+ norm_err_grd2 = 0.0;
+ for (int j = 0; j < dim; ++j)
+ norm_err_grd2 +=
+ sqr(((*grdu)(*worldCoordsAtQP))[j] - grdUhAtQPs[iq][j]);
- val += q->getWeight(iq) * norm_err_grd2;
+ val += q->getWeight(iq) * norm_err_grd2;
- if (relErr) {
- norm_grd2 = 0.0;
- for (int j = 0; j < dim; ++j)
- norm_grd2 += sqr(((*grdu)(*worldCoordsAtQP))[j]);
+ if (relErr) {
+ norm_grd2 = 0.0;
+ for (int j = 0; j < dim; ++j)
+ norm_grd2 += sqr(((*grdu)(*worldCoordsAtQP))[j]);
- val_nrm += q->getWeight(iq) * norm_grd2;
+ val_nrm += q->getWeight(iq) * norm_grd2;
+ }
}
- }
- double err = det * val;
+ double err = det * val;
- if (relErr)
- nrmGrdU += fac * det * val_nrm;
+ if (relErr)
+ nrmGrdU += fac * det * val_nrm;
- return err;
-}
-
-double
-CFE_NormAndErrorFcts::Norm_IntNoBound(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::Norm_IntNoBound()");
-
- int dim = elLS->getDim();
- Mesh *mesh = elLS->getMesh();
- double nrm = 0.0;
- int elStatus;
-
- // ===== Get quadratures. =====
- if (!q) {
- q = Quadrature::provideQuadrature(dim, deg);
+ return err;
}
- elNorm->setQuadrature(q);
+ double
+ CFE_NormAndErrorFcts::Norm_IntNoBound(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::Norm_IntNoBound()");
+
+ int dim = elLS->getDim();
+ Mesh *mesh = elLS->getMesh();
+ double nrm = 0.0;
+ int elStatus;
+
+ // ===== Get quadratures. =====
+ if (!q) {
+ q = Quadrature::provideQuadrature(dim, deg);
+ }
+ elNorm->setQuadrature(q);
- // ===== Traverse mesh and calculate integral on each element. =====
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
+ // ===== Traverse mesh and calculate integral on each element. =====
+ TraverseStack stack;
- while(elInfo) {
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
- // Check whether current element is cut by the zero level set.
- elStatus = elLS->createElementLevelSet(elInfo);
+ while(elInfo) {
- if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
+ // Check whether current element is cut by the zero level set.
+ elStatus = elLS->createElementLevelSet(elInfo);
- // -------------------------------------------------------------------
- // Element is in the domain with negative level set function values.
- // -------------------------------------------------------------------
+ if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ // -------------------------------------------------------------------
+ // Element is in the domain with negative level set function values.
+ // -------------------------------------------------------------------
- nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
- }
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
- elInfo = stack.traverseNext(elInfo);
+ nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
+ }
- } // end of: mesh traverse
+ elInfo = stack.traverseNext(elInfo);
- return nrm;
-}
+ } // end of: mesh traverse
-double
-CFE_NormAndErrorFcts::Norm_IntBound(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::Norm_IntBound()");
-
- int dim = elLS->getDim();
- Mesh *mesh = elLS->getMesh();
- double nrm = 0.0;
- double el_norm;
- VectorOfFixVecs<DimVec<double> > *intersecPts;
- int numIntersecPts;
- SubPolytope *subPolytope;
- ScalableQuadrature *scalQuad;
- int nScalQPts;
- int elStatus;
-
- // ===== Get quadratures. =====
- if (!q) {
- q = Quadrature::provideQuadrature(dim, deg);
+ return nrm;
}
- scalQuad = NEW ScalableQuadrature(q);
- nScalQPts = scalQuad->getNumPoints();
+ double
+ CFE_NormAndErrorFcts::Norm_IntBound(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::Norm_IntBound()");
+
+ int dim = elLS->getDim();
+ Mesh *mesh = elLS->getMesh();
+ double nrm = 0.0;
+ double el_norm;
+ VectorOfFixVecs<DimVec<double> > *intersecPts;
+ int numIntersecPts;
+ SubPolytope *subPolytope;
+ ScalableQuadrature *scalQuad;
+ int nScalQPts;
+ int elStatus;
+
+ // ===== Get quadratures. =====
+ if (!q) {
+ q = Quadrature::provideQuadrature(dim, deg);
+ }
+ scalQuad = NEW ScalableQuadrature(q);
+ nScalQPts = scalQuad->getNumPoints();
- // ===== Traverse mesh and calculate integral on each element. =====
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
+ // ===== Traverse mesh and calculate integral on each element. =====
+ TraverseStack stack;
- while(elInfo) {
- el_norm = 0.0;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
- // Check whether current element is cut by the zero level set.
- elStatus = elLS->createElementLevelSet(elInfo);
+ while(elInfo) {
+ el_norm = 0.0;
- if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+ // Check whether current element is cut by the zero level set.
+ elStatus = elLS->createElementLevelSet(elInfo);
- // -------------------------------------------------------------------
- // Element is cut by the zero level set.
- // -------------------------------------------------------------------
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
- // Calculate norm on subpolyope.
- intersecPts = elLS->getElIntersecPoints();
- numIntersecPts = elLS->getNumElIntersecPoints();
+ // -------------------------------------------------------------------
+ // Element is cut by the zero level set.
+ // -------------------------------------------------------------------
- // -----------------------------------------------------------------
- // Subelement may be inside the domain with negative level set
- // function value as well as inside the domain with positive
- // function value.
- //
- // Whether a subelement is in the domain with negative or positive
- // level set function values is checked by the level set function
- // value of the first vertex of the subelement. (The subelements
- // 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
- // function value of this vertex really is unequal to zero.)
+ // Calculate norm on subpolyope.
+ intersecPts = elLS->getElIntersecPoints();
+ numIntersecPts = elLS->getNumElIntersecPoints();
- subPolytope = NEW SubPolytope(elInfo,
- intersecPts,
- numIntersecPts,
- 0);
+ // -----------------------------------------------------------------
+ // Subelement may be inside the domain with negative level set
+ // function value as well as inside the domain with positive
+ // function value.
+ //
+ // Whether a subelement is in the domain with negative or positive
+ // level set function values is checked by the level set function
+ // value of the first vertex of the subelement. (The subelements
+ // 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
+ // function value of this vertex really is unequal to zero.)
+
+ subPolytope = NEW SubPolytope(elInfo,
+ intersecPts,
+ numIntersecPts,
+ 0);
- elLS->setLevelSetDomain(
- elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
+ elLS->setLevelSetDomain(
+ elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
- el_norm = calcSubPolNorm(elInfo,
- subPolytope,
- elNorm,
- scalQuad);
+ el_norm = calcSubPolNorm(elInfo,
+ subPolytope,
+ elNorm,
+ scalQuad);
- // Calculate integral on the other element part.
- if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR)
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
- else
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
+ // Calculate integral on the other element part.
+ if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR)
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ else
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
- elNorm->setQuadrature(q);
- el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
+ elNorm->setQuadrature(q);
+ el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
- el_norm -= calcSubPolNorm(elInfo,
- subPolytope,
- elNorm,
- scalQuad,
- -1.0);
+ el_norm -= calcSubPolNorm(elInfo,
+ subPolytope,
+ elNorm,
+ scalQuad,
+ -1.0);
- nrm += el_norm;
+ nrm += el_norm;
- // Free data.
- DELETE subPolytope;
- }
- else if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
-
- // -------------------------------------------------------------------
- // Element is in the domain with negative level set function values.
- // -------------------------------------------------------------------
+ // Free data.
+ DELETE subPolytope;
+ }
+ else if (elStatus == ElementLevelSet::LEVEL_SET_INTERIOR) {
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ // -------------------------------------------------------------------
+ // Element is in the domain with negative level set function values.
+ // -------------------------------------------------------------------
- elNorm->setQuadrature(q);
- nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
- }
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
- elInfo = stack.traverseNext(elInfo);
+ elNorm->setQuadrature(q);
+ nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
+ }
- } // end of: mesh traverse
+ elInfo = stack.traverseNext(elInfo);
- DELETE scalQuad;
+ } // end of: mesh traverse
- return nrm;
-}
+ DELETE scalQuad;
-double
-CFE_NormAndErrorFcts::Norm_Int(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::Norm_Int()");
-
- int dim = elLS->getDim();
- Mesh *mesh = elLS->getMesh();
- double nrm = 0.0;
- double el_norm;
- VectorOfFixVecs<DimVec<double> > *intersecPts;
- int numIntersecPts;
- int vertex_interior;
- SubPolytope *subPolytope;
- ScalableQuadrature *scalQuad;
- int nScalQPts;
- int elStatus;
-
- // ===== Get quadratures. =====
- if (!q) {
- q = Quadrature::provideQuadrature(dim, deg);
+ return nrm;
}
- scalQuad = NEW ScalableQuadrature(q);
- nScalQPts = scalQuad->getNumPoints();
+ double
+ CFE_NormAndErrorFcts::Norm_Int(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::Norm_Int()");
+
+ int dim = elLS->getDim();
+ Mesh *mesh = elLS->getMesh();
+ double nrm = 0.0;
+ double el_norm;
+ VectorOfFixVecs<DimVec<double> > *intersecPts;
+ int numIntersecPts;
+ int vertex_interior;
+ SubPolytope *subPolytope;
+ ScalableQuadrature *scalQuad;
+ int nScalQPts;
+ int elStatus;
+
+ // ===== Get quadratures. =====
+ if (!q) {
+ q = Quadrature::provideQuadrature(dim, deg);
+ }
+ scalQuad = NEW ScalableQuadrature(q);
+ nScalQPts = scalQuad->getNumPoints();
- // ===== Traverse mesh and calculate integral on each element. =====
- TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
+ // ===== Traverse mesh and calculate integral on each element. =====
+ TraverseStack stack;
- while(elInfo) {
- el_norm = 0.0;
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
- // Check whether current element is cut by the zero level set.
- elStatus = elLS->createElementLevelSet(elInfo);
+ while(elInfo) {
+ el_norm = 0.0;
- if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+ // Check whether current element is cut by the zero level set.
+ elStatus = elLS->createElementLevelSet(elInfo);
- // -------------------------------------------------------------------
- // Element is cut by the zero level set.
- // -------------------------------------------------------------------
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
- // Create subelements.
- intersecPts = elLS->getElIntersecPoints();
- numIntersecPts = elLS->getNumElIntersecPoints();
+ // -------------------------------------------------------------------
+ // Element is cut by the zero level set.
+ // -------------------------------------------------------------------
- if (dim == 1 || (dim == 3 && numIntersecPts == 4)) {
+ // Create subelements.
+ intersecPts = elLS->getElIntersecPoints();
+ numIntersecPts = elLS->getNumElIntersecPoints();
- // -----------------------------------------------------------------
- // Subelement(s) are inside the domain with negative level set
- // function value.
+ if (dim == 1 || (dim == 3 && numIntersecPts == 4)) {
+
+ // -----------------------------------------------------------------
+ // Subelement(s) are inside the domain with negative level set
+ // function value.
- // Get vertex with negative level set function value.
- for (int i=0; i<=dim; ++i) {
- if (elLS->getElVertLevelSetVec(i) < 0) {
- vertex_interior = i;
- break;
+ // Get vertex with negative level set function value.
+ for (int i=0; i<=dim; ++i) {
+ if (elLS->getElVertLevelSetVec(i) < 0) {
+ vertex_interior = i;
+ break;
+ }
}
+
+ subPolytope = NEW SubPolytope(elInfo,
+ intersecPts,
+ numIntersecPts,
+ vertex_interior);
+
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ }
+ else {
+
+ // -----------------------------------------------------------------
+ // Subelement may be inside the domain with negative level set
+ // function value as well as inside the domain with positive
+ // function value.
+ //
+ // Whether a subelement is in the domain with negative or positive
+ // level set function values is checked by the level set function
+ // value of the first vertex of the subelement. (The subelements
+ // 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
+ // function value of this vertex really is unequal to zero.)
+
+ subPolytope = NEW SubPolytope(elInfo,
+ intersecPts,
+ numIntersecPts,
+ 0);
+
+ elLS->setLevelSetDomain(
+ elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
}
- subPolytope = NEW SubPolytope(elInfo,
- intersecPts,
- numIntersecPts,
- vertex_interior);
+ // Calculate norm on subpolytope.
+ if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_INTERIOR)
+ el_norm = calcSubPolNorm(elInfo,
+ subPolytope,
+ elNorm,
+ scalQuad);
+ else
+ el_norm = calcSubPolNorm(elInfo,
+ subPolytope,
+ elNorm,
+ scalQuad,
+ -1.0);
+
+ // -------------------------------------------------------------------
+ // 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 (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR) {
+
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+
+ elNorm->setQuadrature(q);
+ el_norm *= -1.0;
+ el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
+ }
+
+ // 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);
+
+ elNorm->setQuadrature(q);
+ el_norm = elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
}
- else {
+
+ nrm += el_norm;
+ elInfo = stack.traverseNext(elInfo);
+
+ } // end of: mesh traverse
+
+ DELETE scalQuad;
+
+ return nrm;
+ }
+
+ double
+ CFE_NormAndErrorFcts::Norm_Bound(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::Norm_Bound()");
+
+ int dim = elLS->getDim();
+ Mesh *mesh = elLS->getMesh();
+ double nrm = 0.0;
+ double el_norm;
+ VectorOfFixVecs<DimVec<double> > *intersecPts;
+ int numIntersecPts;
+ SubPolytope *subPolytope;
+ ScalableQuadrature *scalQuad;
+ int nScalQPts;
+ int elStatus;
+
+ // ===== Get quadratures. =====
+ if (!q) {
+ q = Quadrature::provideQuadrature(dim, deg);
+ }
+ scalQuad = NEW ScalableQuadrature(q);
+ nScalQPts = scalQuad->getNumPoints();
+
+
+ // ===== Traverse mesh and calculate integral on each element. =====
+ TraverseStack stack;
+
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
+
+ while(elInfo) {
+ el_norm = 0.0;
+
+ // Check whether current element is cut by the zero level set.
+ elStatus = elLS->createElementLevelSet(elInfo);
+
+ if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
+
+ // -------------------------------------------------------------------
+ // Element is cut by the zero level set.
+ // -------------------------------------------------------------------
+
+ // Calculate norm on subpolyope.
+ intersecPts = elLS->getElIntersecPoints();
+ numIntersecPts = elLS->getNumElIntersecPoints();
// -----------------------------------------------------------------
// Subelement may be inside the domain with negative level set
@@ -457,692 +580,571 @@ CFE_NormAndErrorFcts::Norm_Int(ElementNorm *elNorm,
0);
elLS->setLevelSetDomain(
- elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
- }
+ elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
- // Calculate norm on subpolytope.
- if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_INTERIOR)
el_norm = calcSubPolNorm(elInfo,
subPolytope,
elNorm,
scalQuad);
- else
- el_norm = calcSubPolNorm(elInfo,
- subPolytope,
- elNorm,
- scalQuad,
- -1.0);
- // -------------------------------------------------------------------
- // 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 (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR) {
-
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ // Calculate integral on the other element part.
+ if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR)
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ else
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
elNorm->setQuadrature(q);
- el_norm *= -1.0;
el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
- }
-
- // 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);
-
- elNorm->setQuadrature(q);
- el_norm = elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
- }
+ el_norm -= calcSubPolNorm(elInfo,
+ subPolytope,
+ elNorm,
+ scalQuad,
+ -1.0);
- nrm += el_norm;
- elInfo = stack.traverseNext(elInfo);
+ nrm += el_norm;
+
+ // Free data.
+ DELETE subPolytope;
+ }
- } // end of: mesh traverse
+ elInfo = stack.traverseNext(elInfo);
- DELETE scalQuad;
+ } // end of: mesh traverse
- return nrm;
-}
+ DELETE scalQuad;
-double
-CFE_NormAndErrorFcts::Norm_Bound(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::Norm_Bound()");
-
- int dim = elLS->getDim();
- Mesh *mesh = elLS->getMesh();
- double nrm = 0.0;
- double el_norm;
- VectorOfFixVecs<DimVec<double> > *intersecPts;
- int numIntersecPts;
- SubPolytope *subPolytope;
- ScalableQuadrature *scalQuad;
- int nScalQPts;
- int elStatus;
-
- // ===== Get quadratures. =====
- if (!q) {
- q = Quadrature::provideQuadrature(dim, deg);
+ return nrm;
}
- scalQuad = NEW ScalableQuadrature(q);
- nScalQPts = scalQuad->getNumPoints();
-
-
- // ===== Traverse mesh and calculate integral on each element. =====
- TraverseStack stack;
-
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
-
- while(elInfo) {
- el_norm = 0.0;
-
- // Check whether current element is cut by the zero level set.
- elStatus = elLS->createElementLevelSet(elInfo);
- if (elStatus == ElementLevelSet::LEVEL_SET_BOUNDARY) {
-
- // -------------------------------------------------------------------
- // Element is cut by the zero level set.
- // -------------------------------------------------------------------
-
- // Calculate norm on subpolyope.
- intersecPts = elLS->getElIntersecPoints();
- numIntersecPts = elLS->getNumElIntersecPoints();
-
- // -----------------------------------------------------------------
- // Subelement may be inside the domain with negative level set
- // function value as well as inside the domain with positive
- // function value.
- //
- // Whether a subelement is in the domain with negative or positive
- // level set function values is checked by the level set function
- // value of the first vertex of the subelement. (The subelements
- // 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
- // function value of this vertex really is unequal to zero.)
-
- subPolytope = NEW SubPolytope(elInfo,
- intersecPts,
- numIntersecPts,
- 0);
-
- elLS->setLevelSetDomain(
- elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
-
- el_norm = calcSubPolNorm(elInfo,
- subPolytope,
- elNorm,
- scalQuad);
-
- // Calculate integral on the other element part.
- if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR)
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
- else
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
-
- elNorm->setQuadrature(q);
- el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
-
- el_norm -= calcSubPolNorm(elInfo,
- subPolytope,
- elNorm,
- scalQuad,
- -1.0);
-
- nrm += el_norm;
-
- // Free data.
- DELETE subPolytope;
+ double
+ CFE_NormAndErrorFcts::Norm_Complete(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::Norm_Complete()");
+
+ int dim = elLS->getDim();
+ Mesh *mesh = elLS->getMesh();
+ double nrm = 0.0;
+ double el_norm;
+ VectorOfFixVecs<DimVec<double> > *intersecPts;
+ int numIntersecPts;
+ SubPolytope *subPolytope;
+ ScalableQuadrature *scalQuad;
+ int nScalQPts;
+ int elStatus;
+
+ // ===== Get quadratures. =====
+ if (!q) {
+ q = Quadrature::provideQuadrature(dim, deg);
}
-
- elInfo = stack.traverseNext(elInfo);
-
- } // end of: mesh traverse
-
- DELETE scalQuad;
-
- return nrm;
-}
-
-double
-CFE_NormAndErrorFcts::Norm_Complete(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::Norm_Complete()");
-
- int dim = elLS->getDim();
- Mesh *mesh = elLS->getMesh();
- double nrm = 0.0;
- double el_norm;
- VectorOfFixVecs<DimVec<double> > *intersecPts;
- int numIntersecPts;
- SubPolytope *subPolytope;
- ScalableQuadrature *scalQuad;
- int nScalQPts;
- int elStatus;
-
- // ===== Get quadratures. =====
- if (!q) {
- q = Quadrature::provideQuadrature(dim, deg);
- }
- scalQuad = NEW ScalableQuadrature(q);
- nScalQPts = scalQuad->getNumPoints();
+ scalQuad = NEW ScalableQuadrature(q);
+ nScalQPts = scalQuad->getNumPoints();
- // ===== Traverse mesh and calculate integral on each element. =====
- TraverseStack stack;
+ // ===== Traverse mesh and calculate integral on each element. =====
+ TraverseStack stack;
- ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
+ ElInfo *elInfo = stack.traverseFirst(mesh, -1, fillFlag);
- while(elInfo) {
- el_norm = 0.0;
+ while(elInfo) {
+ el_norm = 0.0;
- // Check whether current element is cut by the zero level set.
- elStatus = elLS->createElementLevelSet(elInfo);
+ // Check whether current element is cut by the zero level set.
+ elStatus = elLS->createElementLevelSet(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.
+ // -------------------------------------------------------------------
- // Calculate norm on subpolyope.
- intersecPts = elLS->getElIntersecPoints();
- numIntersecPts = elLS->getNumElIntersecPoints();
+ // Calculate norm on subpolyope.
+ intersecPts = elLS->getElIntersecPoints();
+ numIntersecPts = elLS->getNumElIntersecPoints();
- // -----------------------------------------------------------------
- // Subelement may be inside the domain with negative level set
- // function value as well as inside the domain with positive
- // function value.
- //
- // Whether a subelement is in the domain with negative or positive
- // level set function values is checked by the level set function
- // value of the first vertex of the subelement. (The subelements
- // 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
- // function value of this vertex really is unequal to zero.)
+ // -----------------------------------------------------------------
+ // Subelement may be inside the domain with negative level set
+ // function value as well as inside the domain with positive
+ // function value.
+ //
+ // Whether a subelement is in the domain with negative or positive
+ // level set function values is checked by the level set function
+ // value of the first vertex of the subelement. (The subelements
+ // 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
+ // function value of this vertex really is unequal to zero.)
- subPolytope = NEW SubPolytope(elInfo,
- intersecPts,
- numIntersecPts,
- 0);
+ subPolytope = NEW SubPolytope(elInfo,
+ intersecPts,
+ numIntersecPts,
+ 0);
- elLS->setLevelSetDomain(
- elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
+ elLS->setLevelSetDomain(
+ elLS->getVertexPos(subPolytope->getSubElement(0)->getLambda(0)));
- el_norm = calcSubPolNorm(elInfo,
- subPolytope,
- elNorm,
- scalQuad);
+ el_norm = calcSubPolNorm(elInfo,
+ subPolytope,
+ elNorm,
+ scalQuad);
- // Calculate integral on the other element part.
- if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR)
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
- else
- elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
+ // Calculate integral on the other element part.
+ if (elLS->getLevelSetDomain() == ElementLevelSet::LEVEL_SET_EXTERIOR)
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_INTERIOR);
+ else
+ elLS->setLevelSetDomain(ElementLevelSet::LEVEL_SET_EXTERIOR);
- elNorm->setQuadrature(q);
- el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
+ elNorm->setQuadrature(q);
+ el_norm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
- el_norm -= calcSubPolNorm(elInfo,
- subPolytope,
- elNorm,
- scalQuad,
- -1.0);
+ el_norm -= calcSubPolNorm(elInfo,
+ subPolytope,
+ elNorm,
+ scalQuad,
+ -1.0);
- nrm += el_norm;
+ nrm += el_norm;
- // Free data.
- DELETE subPolytope;
- }
- else {
-
- // -------------------------------------------------------------------
- // Element is either completely in the domain with negative
- // or positive level set function values.
- // -------------------------------------------------------------------
+ // Free data.
+ DELETE subPolytope;
+ }
+ else {
- elLS->setLevelSetDomain(elStatus);
+ // -------------------------------------------------------------------
+ // Element is either completely in the domain with negative
+ // or positive level set function values.
+ // -------------------------------------------------------------------
- elNorm->setQuadrature(q);
- nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
- }
+ elLS->setLevelSetDomain(elStatus);
- elInfo = stack.traverseNext(elInfo);
+ elNorm->setQuadrature(q);
+ nrm += elNorm->calcElNorm(elInfo, fabs(elInfo->getDet()));
+ }
- } // end of: mesh traverse
+ elInfo = stack.traverseNext(elInfo);
- DELETE scalQuad;
+ } // end of: mesh traverse
- return nrm;
-}
+ DELETE scalQuad;
-double
-CFE_NormAndErrorFcts::L1Norm_Analyt(
- AbstractFunction<double, WorldVector<double> > *f,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::L1Norm_Analyt");
-
- ElementL1Norm_Analyt *elNorm = NEW ElementL1Norm_Analyt(q, f);
- int dim = elLS->getDim();
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET;
-
- double nrm;
- switch(domainFlag) {
- case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ return nrm;
}
- DELETE elNorm;
+ double
+ CFE_NormAndErrorFcts::L1Norm_Analyt(
+ AbstractFunction<double, WorldVector<double> > *f,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::L1Norm_Analyt");
+
+ ElementL1Norm_Analyt *elNorm = NEW ElementL1Norm_Analyt(q, f);
+ int dim = elLS->getDim();
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET;
+
+ double nrm = 0.0;
+ switch(domainFlag) {
+ case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- return nrm;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::L2NormSquare_Analyt(
- AbstractFunction<double, WorldVector<double> > *f,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::L2NormSquare_Analyt");
-
- ElementL2Norm_Analyt *elNorm = NEW ElementL2Norm_Analyt(q, f);
- int dim = elLS->getDim();
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET;
-
- double nrm;
- switch(domainFlag) {
- case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ return nrm;
}
- DELETE elNorm;
+ double
+ CFE_NormAndErrorFcts::L2NormSquare_Analyt(
+ AbstractFunction<double, WorldVector<double> > *f,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::L2NormSquare_Analyt");
+
+ ElementL2Norm_Analyt *elNorm = NEW ElementL2Norm_Analyt(q, f);
+ int dim = elLS->getDim();
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET;
+
+ double nrm = 0.0;
+ switch(domainFlag) {
+ case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- return nrm;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::L2Norm_Analyt(
- AbstractFunction<double, WorldVector<double> > *f,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- return sqrt(L2NormSquare_Analyt(f, elLS, domainFlag, deg, q));
-}
+ return nrm;
+ }
-double
-CFE_NormAndErrorFcts::H1NormSquare_Analyt(
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::H1NormSquare_Analyt");
-
- int dim = elLS->getDim();
- ElementH1Norm_Analyt *elNorm = NEW ElementH1Norm_Analyt(q, grd, dim);
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET;
-
- double nrm;
- switch(domainFlag) {
- case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ double
+ CFE_NormAndErrorFcts::L2Norm_Analyt(
+ AbstractFunction<double, WorldVector<double> > *f,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ return sqrt(L2NormSquare_Analyt(f, elLS, domainFlag, deg, q));
}
- DELETE elNorm;
+ double
+ CFE_NormAndErrorFcts::H1NormSquare_Analyt(
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::H1NormSquare_Analyt");
+
+ int dim = elLS->getDim();
+ ElementH1Norm_Analyt *elNorm = NEW ElementH1Norm_Analyt(q, grd, dim);
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET;
+
+ double nrm = 0.0;
+ switch(domainFlag) {
+ case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- return nrm;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::H1Norm_Analyt(
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- return sqrt(H1NormSquare_Analyt(grd, elLS, domainFlag, deg, q));
-}
+ return nrm;
+ }
-double
-CFE_NormAndErrorFcts::L1Norm_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::L1Norm_DOF");
-
- ElementL1Norm_DOF *elNorm = NEW ElementL1Norm_DOF(q, dof);
- int dim = elLS->getDim();
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET;
-
- double nrm;
- switch(domainFlag) {
- case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ double
+ CFE_NormAndErrorFcts::H1Norm_Analyt(
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ return sqrt(H1NormSquare_Analyt(grd, elLS, domainFlag, deg, q));
}
- DELETE elNorm;
+ double
+ CFE_NormAndErrorFcts::L1Norm_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::L1Norm_DOF");
+
+ ElementL1Norm_DOF *elNorm = NEW ElementL1Norm_DOF(q, dof);
+ int dim = elLS->getDim();
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET;
+
+ double nrm = 0.0;
+ switch(domainFlag) {
+ case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- return nrm;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::L2NormSquare_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::L2NormSquare_DOF");
-
- ElementL2Norm_DOF *elNorm = NEW ElementL2Norm_DOF(q, dof);
- int dim = elLS->getDim();
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET;
-
- double nrm;
- switch(domainFlag) {
- case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ return nrm;
}
- DELETE elNorm;
+ double
+ CFE_NormAndErrorFcts::L2NormSquare_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::L2NormSquare_DOF");
+
+ ElementL2Norm_DOF *elNorm = NEW ElementL2Norm_DOF(q, dof);
+ int dim = elLS->getDim();
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET;
+
+ double nrm = 0.0;
+ switch(domainFlag) {
+ case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- return nrm;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::L2Norm_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- return sqrt(L2NormSquare_DOF(dof, elLS, domainFlag, deg, q));
-}
+ return nrm;
+ }
-double
-CFE_NormAndErrorFcts::H1NormSquare_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::H1NormSquare_DOF");
-
- int dim = elLS->getDim();
- ElementH1Norm_DOF *elNorm = NEW ElementH1Norm_DOF(q, dof, dim);
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA;
-
- double nrm;
- switch(domainFlag) {
- case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ double
+ CFE_NormAndErrorFcts::L2Norm_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ return sqrt(L2NormSquare_DOF(dof, elLS, domainFlag, deg, q));
}
- DELETE elNorm;
+ double
+ CFE_NormAndErrorFcts::H1NormSquare_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::H1NormSquare_DOF");
+
+ int dim = elLS->getDim();
+ ElementH1Norm_DOF *elNorm = NEW ElementH1Norm_DOF(q, dof, dim);
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA;
+
+ double nrm = 0.0;
+ switch(domainFlag) {
+ case -3: nrm = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: nrm = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: nrm = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: nrm = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: nrm = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- return nrm;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::H1Norm_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg,
- Quadrature *q)
-{
- return sqrt(H1NormSquare_DOF(dof, elLS, domainFlag, deg, q));
-}
+ return nrm;
+ }
-double
-CFE_NormAndErrorFcts::L2Err(
- AbstractFunction<double, WorldVector<double> > *u,
- DOFVector<double> *uh,
- ElementLevelSet *elLS,
- int domainFlag,
- int relErr,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::L2Err()");
-
- ElementL2Err *elNorm = NEW ElementL2Err(q, u, uh, relErr);
- int dim = elLS->getDim();
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET;
-
- double err;
- switch(domainFlag) {
- case -3: err = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: err = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: err = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: err = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: err = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ double
+ CFE_NormAndErrorFcts::H1Norm_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg,
+ Quadrature *q)
+ {
+ return sqrt(H1NormSquare_DOF(dof, elLS, domainFlag, deg, q));
}
- L2_err_abs = sqrt(err);
- L2_u_norm = sqrt(elNorm->getNormU());
+ double
+ CFE_NormAndErrorFcts::L2Err(
+ AbstractFunction<double, WorldVector<double> > *u,
+ DOFVector<double> *uh,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int relErr,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::L2Err()");
+
+ ElementL2Err *elNorm = NEW ElementL2Err(q, u, uh, relErr);
+ int dim = elLS->getDim();
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET;
+
+ double err = 0.0;
+ switch(domainFlag) {
+ case -3: err = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: err = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: err = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: err = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: err = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- if (relErr)
- err = L2_err_abs / (L2_u_norm + 1.e-15);
- else
- err = L2_err_abs;
+ L2_err_abs = sqrt(err);
+ L2_u_norm = sqrt(elNorm->getNormU());
- DELETE elNorm;
+ if (relErr)
+ err = L2_err_abs / (L2_u_norm + 1.e-15);
+ else
+ err = L2_err_abs;
- return err;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::H1Err(
- AbstractFunction<WorldVector<double>, WorldVector<double> > *u,
- DOFVector<double> *uh,
- ElementLevelSet *elLS,
- int domainFlag,
- int relErr,
- int deg,
- Quadrature *q)
-{
- FUNCNAME("CFE_NormAndErrorFcts::H1Err()");
-
- int dim = elLS->getDim();
- ElementH1Err *elNorm = NEW ElementH1Err(q, u, uh, relErr, dim);
-
- TEST_EXIT(dim == Global::getGeo(WORLD))
- ("doesn't work for dimension of problem != dimension of world!\n");
-
- Flag fillFlag = Mesh::CALL_LEAF_EL |
- Mesh::FILL_COORDS |
- Mesh::FILL_DET |
- Mesh::FILL_GRD_LAMBDA;
-
- double err;
- switch(domainFlag) {
- case -3: err = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -2: err = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
- break;
- case -1: err = Norm_Int(elNorm, elLS, fillFlag, deg, q);
- break;
- case 0: err = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
- break;
- case 1: err = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
- break;
- default: ERROR_EXIT("illegal flag !\n");
- break;
+ return err;
}
- H1_err_abs = sqrt(err);
- H1_u_norm = sqrt(elNorm->getNormGrdU());
+ double
+ CFE_NormAndErrorFcts::H1Err(
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *u,
+ DOFVector<double> *uh,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int relErr,
+ int deg,
+ Quadrature *q)
+ {
+ FUNCNAME("CFE_NormAndErrorFcts::H1Err()");
+
+ int dim = elLS->getDim();
+ ElementH1Err *elNorm = NEW ElementH1Err(q, u, uh, relErr, dim);
+
+ TEST_EXIT(dim == Global::getGeo(WORLD))
+ ("doesn't work for dimension of problem != dimension of world!\n");
+
+ Flag fillFlag = Mesh::CALL_LEAF_EL |
+ Mesh::FILL_COORDS |
+ Mesh::FILL_DET |
+ Mesh::FILL_GRD_LAMBDA;
+
+ double err = 0.0;
+ switch(domainFlag) {
+ case -3: err = Norm_IntNoBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -2: err = Norm_IntBound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case -1: err = Norm_Int(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 0: err = Norm_Bound(elNorm, elLS, fillFlag, deg, q);
+ break;
+ case 1: err = Norm_Complete(elNorm, elLS, fillFlag, deg, q);
+ break;
+ default: ERROR_EXIT("illegal flag !\n");
+ break;
+ }
- if (relErr)
- err = H1_err_abs / (H1_u_norm + 1.e-15);
- else
- err = H1_err_abs;
+ H1_err_abs = sqrt(err);
+ H1_u_norm = sqrt(elNorm->getNormGrdU());
- DELETE elNorm;
+ if (relErr)
+ err = H1_err_abs / (H1_u_norm + 1.e-15);
+ else
+ err = H1_err_abs;
- return err;
-}
+ DELETE elNorm;
-double
-CFE_NormAndErrorFcts::calcSubPolNorm(ElInfo *elInfo,
- SubPolytope *subPolytope,
- ElementNorm *elNorm,
- ScalableQuadrature *scalQuad,
- const double &subPolFac)
-{
- double nrm = 0.0;
-
- for (std::vector<SubElInfo *>::iterator it =
- subPolytope->getSubElementsBegin();
- it != subPolytope->getSubElementsEnd();
- it++) {
-
- scalQuad->scaleQuadrature(**it);
- elNorm->setQuadrature(scalQuad);
+ return err;
+ }
+
+ double
+ CFE_NormAndErrorFcts::calcSubPolNorm(ElInfo *elInfo,
+ SubPolytope *subPolytope,
+ ElementNorm *elNorm,
+ ScalableQuadrature *scalQuad,
+ const double &subPolFac)
+ {
+ double nrm = 0.0;
+
+ for (std::vector<SubElInfo *>::iterator it =
+ subPolytope->getSubElementsBegin();
+ it != subPolytope->getSubElementsEnd();
+ it++) {
+
+ scalQuad->scaleQuadrature(**it);
+ elNorm->setQuadrature(scalQuad);
- nrm += elNorm->calcElNorm(elInfo,
- fabs((*it)->getDet()),
- subPolFac);
+ nrm += elNorm->calcElNorm(elInfo,
+ fabs((*it)->getDet()),
+ subPolFac);
+ }
+
+ return nrm;
}
- return nrm;
}
diff --git a/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.h b/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.h
index 4363642422b33c14991b24c4ac0b5db65ec1f118..cd8c616e43d7a15739c9867ca2c27c93558e92bc 100644
--- a/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.h
+++ b/AMDiS/compositeFEM/src/CFE_NormAndErrorFcts.h
@@ -11,564 +11,566 @@
#include "ScalableQuadrature.h"
#include "SubPolytope.h"
-using namespace AMDiS;
-
-class ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementNorm);
-
- /**
- * Constructor.
- */
- ElementNorm(Quadrature *q_)
- : q(q_),
- nQPts(0)
+namespace AMDiS {
+
+ class ElementNorm
{
- if (q)
+ public:
+ MEMORY_MANAGED(ElementNorm);
+
+ /**
+ * Constructor.
+ */
+ ElementNorm(Quadrature *q_)
+ : q(q_),
+ nQPts(0)
+ {
+ if (q)
+ nQPts = q->getNumPoints();
+ };
+
+ /**
+ * Destructor.
+ */
+ virtual ~ElementNorm()
+ {};
+
+ /**
+ * Calculates element norm on elInfo.
+ */
+ virtual double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0) = 0;
+
+ /**
+ * Sets quadrature to q_.
+ */
+ inline void setQuadrature(Quadrature *q_) {
+ q = q_;
nQPts = q->getNumPoints();
+ };
+
+ protected:
+ /**
+ * Quadrature formula.
+ */
+ Quadrature *q;
+
+ /**
+ * Number of quadrature points.
+ */
+ int nQPts;
};
- /**
- * Destructor.
- */
- virtual ~ElementNorm()
- {};
-
- /**
- * Calculates element norm on elInfo.
- */
- virtual double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0) = 0;
-
- /**
- * Sets quadrature to q_.
- */
- inline void setQuadrature(Quadrature *q_) {
- q = q_;
- nQPts = q->getNumPoints();
+ class ElementL1Norm_Analyt : public ElementNorm
+ {
+ public:
+ MEMORY_MANAGED(ElementL1Norm_Analyt);
+
+ /**
+ * Constructor.
+ */
+ ElementL1Norm_Analyt(Quadrature *q_,
+ AbstractFunction<double, WorldVector<double> > *f_)
+ : ElementNorm(q_),
+ f(f_)
+ {};
+
+ /**
+ * Calculates element norm on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ protected:
+ /**
+ * Abstract function for which norm is calculated.
+ */
+ AbstractFunction<double, WorldVector<double> > *f;
};
- protected:
- /**
- * Quadrature formula.
- */
- Quadrature *q;
-
- /**
- * Number of quadrature points.
- */
- int nQPts;
-};
-
-class ElementL1Norm_Analyt : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementL1Norm_Analyt);
-
- /**
- * Constructor.
- */
- ElementL1Norm_Analyt(Quadrature *q_,
- AbstractFunction<double, WorldVector<double> > *f_)
- : ElementNorm(q_),
- f(f_)
- {};
-
- /**
- * Calculates element norm on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- protected:
- /**
- * Abstract function for which norm is calculated.
- */
- AbstractFunction<double, WorldVector<double> > *f;
-};
-
-class ElementL2Norm_Analyt : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementL2Norm_Analyt);
-
- /**
- * Constructor.
- */
- ElementL2Norm_Analyt(Quadrature *q_,
- AbstractFunction<double, WorldVector<double> > *f_)
- : ElementNorm(q_),
- f(f_)
- {};
-
- /**
- * Calculates element norm on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- protected:
- /**
- * Abstract function for which norm is calculated.
- */
- AbstractFunction<double, WorldVector<double> > *f;
-};
-
-class ElementH1Norm_Analyt : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementH1Norm_Analyt);
-
- /**
- * Constructor.
- */
- ElementH1Norm_Analyt(Quadrature *q_,
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grd_,
- int dim_)
- : ElementNorm(q_),
- grd(grd_),
- dim(dim_)
- {};
-
- /**
- * Calculates element norm on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- protected:
- /**
- * Abstract function for which norm is calculated.
- */
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grd;
-
- /**
- * Mesh dimension.
- */
- int dim;
-};
-
-class ElementL1Norm_DOF : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementL1Norm_DOF);
-
- /**
- * Constructor.
- */
- ElementL1Norm_DOF(Quadrature *q_, DOFVector<double> *dofVec_)
- : ElementNorm(q_),
- dofVec(dofVec_)
- {};
-
- /**
- * Calculates element norm on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- protected:
- /**
- * DOF vector for which norm is calculated.
- */
- DOFVector<double> *dofVec;
-};
-
-class ElementL2Norm_DOF : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementL2Norm_DOF);
-
- /**
- * Constructor.
- */
- ElementL2Norm_DOF(Quadrature *q_, DOFVector<double> *dofVec_)
- : ElementNorm(q_),
- dofVec(dofVec_)
- {};
-
- /**
- * Calculates element norm on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- protected:
- /**
- * DOF vector for which norm is calculated.
- */
- DOFVector<double> *dofVec;
-};
-
-class ElementH1Norm_DOF : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementH1Norm_DOF);
-
- /**
- * Constructor.
- */
- ElementH1Norm_DOF(Quadrature *q_, DOFVector<double> *dofVec_, int dim_)
- : ElementNorm(q_),
- dofVec(dofVec_),
- dim(dim_)
- {};
-
- /**
- * Calculates element norm on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- protected:
- /**
- * DOF vector for which norm is calculated.
- */
- DOFVector<double> *dofVec;
-
- /**
- * Mesh dimension.
- */
- int dim;
-};
-
-class ElementL2Err : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementL2Err);
-
- /**
- * Constructor.
- */
- ElementL2Err(Quadrature *q_,
- AbstractFunction<double, WorldVector<double> > *u_,
- DOFVector<double> *uh_,
- int relErr_)
- : ElementNorm(q_),
- u(u_),
- uh(uh_),
- relErr(relErr_),
- nrmU(0.0)
- {};
-
- /**
- * Calculates element error on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- /**
- * Get norm of u.
- */
- inline double getNormU() const {
- return nrmU;
+ class ElementL2Norm_Analyt : public ElementNorm
+ {
+ public:
+ MEMORY_MANAGED(ElementL2Norm_Analyt);
+
+ /**
+ * Constructor.
+ */
+ ElementL2Norm_Analyt(Quadrature *q_,
+ AbstractFunction<double, WorldVector<double> > *f_)
+ : ElementNorm(q_),
+ f(f_)
+ {};
+
+ /**
+ * Calculates element norm on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ protected:
+ /**
+ * Abstract function for which norm is calculated.
+ */
+ AbstractFunction<double, WorldVector<double> > *f;
};
- protected:
- /**
- * Abstract function for which error is calculated.
- */
- AbstractFunction<double, WorldVector<double> > *u;
-
- /**
- * DOF vector for which error is calculated.
- */
- DOFVector<double> *uh;
-
- /**
- * Indicates whether relative (1) or absolute error (0) is calculated.
- */
- int relErr;
-
- /**
- * Norm of u in case relative error is calculated.
- */
- double nrmU;
-};
-
-class ElementH1Err : public ElementNorm
-{
- public:
- MEMORY_MANAGED(ElementH1Err);
-
- /**
- * Constructor.
- */
- ElementH1Err(Quadrature *q_,
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grdu_,
- DOFVector<double> *uh_,
- int relErr_,
- int dim_)
- : ElementNorm(q_),
- grdu(grdu_),
- uh(uh_),
- relErr(relErr_),
- nrmGrdU(0.0),
- dim(dim_)
- {};
-
- /**
- * Calculates element error on elInfo.
- */
- double calcElNorm(ElInfo *elInfo,
- const double &det,
- const double &fac = 1.0);
-
- /**
- * Get norm of grdu.
- */
- inline double getNormGrdU() const {
- return nrmGrdU;
+
+ class ElementH1Norm_Analyt : public ElementNorm
+ {
+ public:
+ MEMORY_MANAGED(ElementH1Norm_Analyt);
+
+ /**
+ * Constructor.
+ */
+ ElementH1Norm_Analyt(Quadrature *q_,
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grd_,
+ int dim_)
+ : ElementNorm(q_),
+ grd(grd_),
+ dim(dim_)
+ {};
+
+ /**
+ * Calculates element norm on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ protected:
+ /**
+ * Abstract function for which norm is calculated.
+ */
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grd;
+
+ /**
+ * Mesh dimension.
+ */
+ int dim;
};
- protected:
- /**
- * Abstract function for which norm is calculated.
- */
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grdu;
-
- /**
- * DOF vector for which error is calculated.
- */
- DOFVector<double> *uh;
-
- /**
- * Indicates whether relative (1) or absolute error (0) is calculated.
- */
- int relErr;
-
- /**
- * Norm of grdu in case relative error is calculated.
- */
- double nrmGrdU;
-
- /**
- * Mesh dimension.
- */
- int dim;
-};
-
-/////////////////////////////////////////////////////////////////////////////
-///// class CFE_NormAndErrorFcts ////////////////////////////////////////////
-/////////////////////////////////////////////////////////////////////////////
-class CFE_NormAndErrorFcts
-{
- public:
- MEMORY_MANAGED(CFE_NormAndErrorFcts);
-
- // ========================================================================
- // Calculation of L1-Norm, L2-Norm and H1-Norm on domain depending
- // on the level set function.
- //
- // Here, the H1-norm is the H1-seminorm, i.e. you get the "common"
- // H1-norm by
- // H1_Norm() = sqrt(L2_Norm^2() + H1_Seminorm^2()) .
- //
- // Parameters for domainFlag:
- // -3: all elements in the domain with negative level set function
- // values which do not intersect the zero level set
- // -2 : all elements which intersect the domain with negative level set
- // function values, in particular all (complete) elements which
- // cut the zero level set
- // -1 : domain with negative level set function values
- // 0 : all elements cut by the zero level set
- // 1 : complete mesh
- // ========================================================================
- static double L1Norm_Analyt(
- AbstractFunction<double, WorldVector<double> > *f,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 1,
- Quadrature* q = NULL);
- static double L2Norm_Analyt(
- AbstractFunction<double, WorldVector<double> > *f,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 2,
- Quadrature* q = NULL);
- static double L2NormSquare_Analyt(
- AbstractFunction<double, WorldVector<double> > *f,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 2,
- Quadrature* q = NULL);
- static double H1Norm_Analyt(
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 0,
- Quadrature* q = NULL);
- static double H1NormSquare_Analyt(
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 0,
- Quadrature* q = NULL);
- static double L1Norm_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 1,
- Quadrature* q = NULL);
- static double L2Norm_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 2,
- Quadrature* q = NULL);
- static double L2NormSquare_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 2,
- Quadrature* q = NULL);
- static double H1Norm_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 0,
- Quadrature* q = NULL);
- static double H1NormSquare_DOF(DOFVector<double> *dof,
- ElementLevelSet *elLS,
- int domainFlag,
- int deg = 0,
- Quadrature* q = NULL);
-
- // ========================================================================
- // Calculation of error between
- // -> true solution (analytically given) and
- // -> FE solution (dof vector).
- //
- // Here, the H1-norm is the H1-seminorm, i.e. you get the "common"
- // H1-norm by
- // H1_Norm() = sqrt(L2_Norm^2() + H1_Seminorm^2()) .
- //
- // Parameter \ref domainFlag:
- // -3: all elements in the domain with negative level set function
- // values which do not intersect the zero level set
- // -2 : all elements which intersect the domain with negative level set
- // function values, in particular all (complete) elements which
- // cut the zero level set
- // -1 : domain with negative level set function values
- // 0 : all elements cut by the zero level set
- // 1 : complete mesh
- // ========================================================================
- static double L2Err(AbstractFunction<double, WorldVector<double> > *u,
- DOFVector<double> *uh,
- ElementLevelSet *elLS,
- int domainFlag,
- int relErr = 0,
- int deg = 2,
- Quadrature *q = NULL);
- static double H1Err(
- AbstractFunction<WorldVector<double>, WorldVector<double> > *grdU,
- DOFVector<double> *uh,
- ElementLevelSet *elLS,
- int domainFlag,
- int relErr = 0,
- int deg = 0,
- Quadrature *q = NULL);
-
- /**
- * Get absolute L2 error.
- * (If relative L2 error is calculated, the absolute L2 error is also
- * calculated and stored in L2_err_abs).
- */
- inline static double getL2ErrAbs()
+ class ElementL1Norm_DOF : public ElementNorm
{
- return L2_err_abs;
- }
-
- /**
- * Get absolute H1 error.
- * (If relative H1 error is calculated, the absolute H1 error is also
- * calculated and stored in H1_err_abs).
- */
- inline static double getH1ErrAbs()
+ public:
+ MEMORY_MANAGED(ElementL1Norm_DOF);
+
+ /**
+ * Constructor.
+ */
+ ElementL1Norm_DOF(Quadrature *q_, DOFVector<double> *dofVec_)
+ : ElementNorm(q_),
+ dofVec(dofVec_)
+ {};
+
+ /**
+ * Calculates element norm on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ protected:
+ /**
+ * DOF vector for which norm is calculated.
+ */
+ DOFVector<double> *dofVec;
+ };
+
+ class ElementL2Norm_DOF : public ElementNorm
{
- return H1_err_abs;
- }
-
- /**
- * Get L2 norm of solution u.
- * (If relative L2 error is calculated, the L2 norm of the solution u is
- * also calculated and stored in L2_u_norm).
- */
- inline static double getL2_U_Norm()
+ public:
+ MEMORY_MANAGED(ElementL2Norm_DOF);
+
+ /**
+ * Constructor.
+ */
+ ElementL2Norm_DOF(Quadrature *q_, DOFVector<double> *dofVec_)
+ : ElementNorm(q_),
+ dofVec(dofVec_)
+ {};
+
+ /**
+ * Calculates element norm on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ protected:
+ /**
+ * DOF vector for which norm is calculated.
+ */
+ DOFVector<double> *dofVec;
+ };
+
+ class ElementH1Norm_DOF : public ElementNorm
{
- return L2_u_norm;
- }
-
- /**
- * Get H1 norm of solution u.
- * (If relative H1 error is calculated, the H1 norm of the solution u is
- * also calculated and stored in H1_u_norm).
- */
- inline static double getH1_U_Norm()
+ public:
+ MEMORY_MANAGED(ElementH1Norm_DOF);
+
+ /**
+ * Constructor.
+ */
+ ElementH1Norm_DOF(Quadrature *q_, DOFVector<double> *dofVec_, int dim_)
+ : ElementNorm(q_),
+ dofVec(dofVec_),
+ dim(dim_)
+ {};
+
+ /**
+ * Calculates element norm on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ protected:
+ /**
+ * DOF vector for which norm is calculated.
+ */
+ DOFVector<double> *dofVec;
+
+ /**
+ * Mesh dimension.
+ */
+ int dim;
+ };
+
+ class ElementL2Err : public ElementNorm
{
- return H1_u_norm;
- }
+ public:
+ MEMORY_MANAGED(ElementL2Err);
+
+ /**
+ * Constructor.
+ */
+ ElementL2Err(Quadrature *q_,
+ AbstractFunction<double, WorldVector<double> > *u_,
+ DOFVector<double> *uh_,
+ int relErr_)
+ : ElementNorm(q_),
+ u(u_),
+ uh(uh_),
+ relErr(relErr_),
+ nrmU(0.0)
+ {};
+
+ /**
+ * Calculates element error on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ /**
+ * Get norm of u.
+ */
+ inline double getNormU() const {
+ return nrmU;
+ };
+ protected:
+ /**
+ * Abstract function for which error is calculated.
+ */
+ AbstractFunction<double, WorldVector<double> > *u;
+
+ /**
+ * DOF vector for which error is calculated.
+ */
+ DOFVector<double> *uh;
+
+ /**
+ * Indicates whether relative (1) or absolute error (0) is calculated.
+ */
+ int relErr;
+
+ /**
+ * Norm of u in case relative error is calculated.
+ */
+ double nrmU;
+ };
-protected:
- static double Norm_IntNoBound(ElementNorm *elNorm,
+ class ElementH1Err : public ElementNorm
+ {
+ public:
+ MEMORY_MANAGED(ElementH1Err);
+
+ /**
+ * Constructor.
+ */
+ ElementH1Err(Quadrature *q_,
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grdu_,
+ DOFVector<double> *uh_,
+ int relErr_,
+ int dim_)
+ : ElementNorm(q_),
+ grdu(grdu_),
+ uh(uh_),
+ relErr(relErr_),
+ nrmGrdU(0.0),
+ dim(dim_)
+ {};
+
+ /**
+ * Calculates element error on elInfo.
+ */
+ double calcElNorm(ElInfo *elInfo,
+ const double &det,
+ const double &fac = 1.0);
+
+ /**
+ * Get norm of grdu.
+ */
+ inline double getNormGrdU() const {
+ return nrmGrdU;
+ };
+
+ protected:
+ /**
+ * Abstract function for which norm is calculated.
+ */
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grdu;
+
+ /**
+ * DOF vector for which error is calculated.
+ */
+ DOFVector<double> *uh;
+
+ /**
+ * Indicates whether relative (1) or absolute error (0) is calculated.
+ */
+ int relErr;
+
+ /**
+ * Norm of grdu in case relative error is calculated.
+ */
+ double nrmGrdU;
+
+ /**
+ * Mesh dimension.
+ */
+ int dim;
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ ///// class CFE_NormAndErrorFcts ////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////
+ class CFE_NormAndErrorFcts
+ {
+ public:
+ MEMORY_MANAGED(CFE_NormAndErrorFcts);
+
+ // ========================================================================
+ // Calculation of L1-Norm, L2-Norm and H1-Norm on domain depending
+ // on the level set function.
+ //
+ // Here, the H1-norm is the H1-seminorm, i.e. you get the "common"
+ // H1-norm by
+ // H1_Norm() = sqrt(L2_Norm^2() + H1_Seminorm^2()) .
+ //
+ // Parameters for domainFlag:
+ // -3: all elements in the domain with negative level set function
+ // values which do not intersect the zero level set
+ // -2 : all elements which intersect the domain with negative level set
+ // function values, in particular all (complete) elements which
+ // cut the zero level set
+ // -1 : domain with negative level set function values
+ // 0 : all elements cut by the zero level set
+ // 1 : complete mesh
+ // ========================================================================
+ static double L1Norm_Analyt(
+ AbstractFunction<double, WorldVector<double> > *f,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 1,
+ Quadrature* q = NULL);
+ static double L2Norm_Analyt(
+ AbstractFunction<double, WorldVector<double> > *f,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 2,
+ Quadrature* q = NULL);
+ static double L2NormSquare_Analyt(
+ AbstractFunction<double, WorldVector<double> > *f,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 2,
+ Quadrature* q = NULL);
+ static double H1Norm_Analyt(
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 0,
+ Quadrature* q = NULL);
+ static double H1NormSquare_Analyt(
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grd,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 0,
+ Quadrature* q = NULL);
+ static double L1Norm_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 1,
+ Quadrature* q = NULL);
+ static double L2Norm_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 2,
+ Quadrature* q = NULL);
+ static double L2NormSquare_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 2,
+ Quadrature* q = NULL);
+ static double H1Norm_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 0,
+ Quadrature* q = NULL);
+ static double H1NormSquare_DOF(DOFVector<double> *dof,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int deg = 0,
+ Quadrature* q = NULL);
+
+ // ========================================================================
+ // Calculation of error between
+ // -> true solution (analytically given) and
+ // -> FE solution (dof vector).
+ //
+ // Here, the H1-norm is the H1-seminorm, i.e. you get the "common"
+ // H1-norm by
+ // H1_Norm() = sqrt(L2_Norm^2() + H1_Seminorm^2()) .
+ //
+ // Parameter \ref domainFlag:
+ // -3: all elements in the domain with negative level set function
+ // values which do not intersect the zero level set
+ // -2 : all elements which intersect the domain with negative level set
+ // function values, in particular all (complete) elements which
+ // cut the zero level set
+ // -1 : domain with negative level set function values
+ // 0 : all elements cut by the zero level set
+ // 1 : complete mesh
+ // ========================================================================
+ static double L2Err(AbstractFunction<double, WorldVector<double> > *u,
+ DOFVector<double> *uh,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int relErr = 0,
+ int deg = 2,
+ Quadrature *q = NULL);
+ static double H1Err(
+ AbstractFunction<WorldVector<double>, WorldVector<double> > *grdU,
+ DOFVector<double> *uh,
+ ElementLevelSet *elLS,
+ int domainFlag,
+ int relErr = 0,
+ int deg = 0,
+ Quadrature *q = NULL);
+
+ /**
+ * Get absolute L2 error.
+ * (If relative L2 error is calculated, the absolute L2 error is also
+ * calculated and stored in L2_err_abs).
+ */
+ inline static double getL2ErrAbs()
+ {
+ return L2_err_abs;
+ }
+
+ /**
+ * Get absolute H1 error.
+ * (If relative H1 error is calculated, the absolute H1 error is also
+ * calculated and stored in H1_err_abs).
+ */
+ inline static double getH1ErrAbs()
+ {
+ return H1_err_abs;
+ }
+
+ /**
+ * Get L2 norm of solution u.
+ * (If relative L2 error is calculated, the L2 norm of the solution u is
+ * also calculated and stored in L2_u_norm).
+ */
+ inline static double getL2_U_Norm()
+ {
+ return L2_u_norm;
+ }
+
+ /**
+ * Get H1 norm of solution u.
+ * (If relative H1 error is calculated, the H1 norm of the solution u is
+ * also calculated and stored in H1_u_norm).
+ */
+ inline static double getH1_U_Norm()
+ {
+ return H1_u_norm;
+ }
+
+ protected:
+ static double Norm_IntNoBound(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature* q);
+ static double Norm_IntBound(ElementNorm *elNorm,
ElementLevelSet *elLS,
Flag fillFlag,
int deg,
Quadrature* q);
- static double Norm_IntBound(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature* q);
- static double Norm_Int(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature* q);
- static double Norm_Bound(ElementNorm *elNorm,
+ static double Norm_Int(ElementNorm *elNorm,
ElementLevelSet *elLS,
Flag fillFlag,
int deg,
Quadrature* q);
- static double Norm_Complete(ElementNorm *elNorm,
- ElementLevelSet *elLS,
- Flag fillFlag,
- int deg,
- Quadrature* q);
-
- /**
- * Calculate norm on subpolytope.
- */
- static double calcSubPolNorm(ElInfo *elInfo,
- SubPolytope *subPolytope,
- ElementNorm *elNorm,
- ScalableQuadrature *scalQuad,
- const double &subPolFac = 1.0);
-
- protected:
- /**
- * Absolute L2 error (last L2 error calculation !).
- */
- static double L2_err_abs;
-
- /**
- * L2 norm of correct solution u (last L2 error calculation !).
- */
- static double L2_u_norm;
-
- /**
- * Absolute H1 error (last H1 error calculation !).
- */
- static double H1_err_abs;
-
- /**
- * H1 norm of correct solution u (last H1 error calculation !).
- */
- static double H1_u_norm;
-};
+ static double Norm_Bound(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature* q);
+ static double Norm_Complete(ElementNorm *elNorm,
+ ElementLevelSet *elLS,
+ Flag fillFlag,
+ int deg,
+ Quadrature* q);
+
+ /**
+ * Calculate norm on subpolytope.
+ */
+ static double calcSubPolNorm(ElInfo *elInfo,
+ SubPolytope *subPolytope,
+ ElementNorm *elNorm,
+ ScalableQuadrature *scalQuad,
+ const double &subPolFac = 1.0);
+
+ protected:
+ /**
+ * Absolute L2 error (last L2 error calculation !).
+ */
+ static double L2_err_abs;
+
+ /**
+ * L2 norm of correct solution u (last L2 error calculation !).
+ */
+ static double L2_u_norm;
+
+ /**
+ * Absolute H1 error (last H1 error calculation !).
+ */
+ static double H1_err_abs;
+
+ /**
+ * H1 norm of correct solution u (last H1 error calculation !).
+ */
+ static double H1_u_norm;
+ };
+
+}
#endif // AMDIS_CFE_NORMANDERRORFCTS_H
diff --git a/AMDiS/compositeFEM/src/ScalableQuadrature.cc b/AMDiS/compositeFEM/src/ScalableQuadrature.cc
index 38e31970ba105b45905cc4dbdd8974a27fe7d7ee..6db04c06df65cd80982a760cd5706e6c495a32b3 100644
--- a/AMDiS/compositeFEM/src/ScalableQuadrature.cc
+++ b/AMDiS/compositeFEM/src/ScalableQuadrature.cc
@@ -1,56 +1,55 @@
#include "ScalableQuadrature.h"
#include "SubElInfo.h"
-ScalableQuadrature::ScalableQuadrature(Quadrature *quadrature)
- : Quadrature(*quadrature)
-{
- /**
- * Change name of quadrature.
- */
- name = "Scalable" + getName();
-
- /**
- * copy quadrature->lambda to oldLambda
- */
- oldLambda = NEW VectorOfFixVecs<DimVec<double> >(*(quadrature->getLambda()));
-}
+namespace AMDiS {
+
+ ScalableQuadrature::ScalableQuadrature(Quadrature *quadrature)
+ : Quadrature(*quadrature)
+ {
+ /**
+ * Change name of quadrature.
+ */
+ name = "Scalable" + getName();
+
+ /**
+ * copy quadrature->lambda to oldLambda
+ */
+ oldLambda = NEW VectorOfFixVecs<DimVec<double> >(*(quadrature->getLambda()));
+ }
-void
-ScalableQuadrature::scaleQuadrature(const SubElInfo& subElInfo)
-{
- //******************************************************************************
- // Manipulates the quadrature points for the assemblage of a subelement.
- // A quadrature point for the assemblage of the subelement is then given in
- // barycentric coordinates with respect to the element containing the
- // subelement.
- // Thus the assemblage of the subelement is done by assembling the element
- // using the quadrature with manipulated quadrature points.
- // Note: The result must be corrected by the determinant of the subelement !!!
- //******************************************************************************
- int iq, i, j;
- double l;
-
- /**
- * If (l_0, l_1, l_2) is the old quadrature point and x_0, x_1, x_2 are the
- * barycentric coordinates of the vertices of the subelement, the new quadrature
- * point (n_0, n_1, n_2) is given as follows:
- * (n_0, n_1, n_2) = l_0 * x_0 + l_1 * x_1 + l_2 * x_2 .
- */
- for (iq = 0; iq < n_points; iq++) {
-
- for (i = 0; i <= dim; i++) {
-
- /**
- * Calculate the i-th component of the iq-th new quadrature point.
- */
- l = 0.0;
-
- for (j = 0; j <= dim; j++) {
- l += getOldLambda(iq, j) * subElInfo.getLambda(j, i);
+ void ScalableQuadrature::scaleQuadrature(const SubElInfo& subElInfo)
+ {
+ //******************************************************************************
+ // Manipulates the quadrature points for the assemblage of a subelement.
+ // A quadrature point for the assemblage of the subelement is then given in
+ // barycentric coordinates with respect to the element containing the
+ // subelement.
+ // Thus the assemblage of the subelement is done by assembling the element
+ // using the quadrature with manipulated quadrature points.
+ // Note: The result must be corrected by the determinant of the subelement !!!
+ //******************************************************************************
+
+ /**
+ * If (l_0, l_1, l_2) is the old quadrature point and x_0, x_1, x_2 are the
+ * barycentric coordinates of the vertices of the subelement, the new quadrature
+ * point (n_0, n_1, n_2) is given as follows:
+ * (n_0, n_1, n_2) = l_0 * x_0 + l_1 * x_1 + l_2 * x_2 .
+ */
+ for (int iq = 0; iq < n_points; iq++) {
+
+ for (int i = 0; i <= dim; i++) {
+
+ /**
+ * Calculate the i-th component of the iq-th new quadrature point.
+ */
+ double l = 0.0;
+
+ for (int j = 0; j <= dim; j++) {
+ l += getOldLambda(iq, j) * subElInfo.getLambda(j, i);
+ }
+
+ (*lambda)[iq][i] = l;
}
-
- (*lambda)[iq][i] = l;
}
}
}
-
diff --git a/AMDiS/compositeFEM/src/ScalableQuadrature.h b/AMDiS/compositeFEM/src/ScalableQuadrature.h
index 4f1fea012f55497d42c86bfae7380b91f1b31c99..7797f6b52272bc88ca64da4827d0b3e46b7fe666 100644
--- a/AMDiS/compositeFEM/src/ScalableQuadrature.h
+++ b/AMDiS/compositeFEM/src/ScalableQuadrature.h
@@ -7,76 +7,78 @@
#include "Quadrature.h"
#include "FixVec.h"
-class SubElInfo;
+namespace AMDiS {
-using namespace AMDiS;
-using namespace std;
+ class SubElInfo;
-// ===============================================================================
-// ===== class ScalableQuadrature ================================================
-// ===============================================================================
-//
-// Class Description:
-// The class \ref ScalableQuadrature holds the functionality for the manipulation
-// of a quadrature formula used for the integration on subelements.
-// If S' is a sublement and S the element containing S', the numerical integration
-// on S' is done by integration on the element S with a manipulated quadrature
-// formula and multiplication of the result with a correction term consisting of
-// the determinants corresponding to S and S'. That means, the quadrature points
-// are manipulated in the following way:
-// The original quadrature formula holds quadrature points which are given in
-// barycentric coordinates with respect to S'. Now we express these quadrature
-// points in barycentric coordinates with respect to S and obtain the manipulated
-// quadrature points we need. Obviously, the corresponding quadrature points
-// in world coordinates coincide. Thus the numerical integration on S with the
-// manipulated quadrature formula gives us the same result as the numerical
-// integration on S' with the original quadrature formula up to the determinant.
-// This method for integration on subelements allows the reuse of the routines for
-// the integration on elements.
-//
-// Main routines:
-// ScalableQuadrature() - Create a new quadrature formula which is a copy of the
-// original quadrature formula and store the original
-// quadrature points in \ref oldLambda.
-// scaleQuadrature() - Manipulate the original quadrature formula for a
-// subelement.
-// ===============================================================================
-class ScalableQuadrature : public Quadrature
-{
-public:
- MEMORY_MANAGED(ScalableQuadrature);
+ // ===============================================================================
+ // ===== class ScalableQuadrature ================================================
+ // ===============================================================================
+ //
+ // Class Description:
+ // The class \ref ScalableQuadrature holds the functionality for the manipulation
+ // of a quadrature formula used for the integration on subelements.
+ // If S' is a sublement and S the element containing S', the numerical integration
+ // on S' is done by integration on the element S with a manipulated quadrature
+ // formula and multiplication of the result with a correction term consisting of
+ // the determinants corresponding to S and S'. That means, the quadrature points
+ // are manipulated in the following way:
+ // The original quadrature formula holds quadrature points which are given in
+ // barycentric coordinates with respect to S'. Now we express these quadrature
+ // points in barycentric coordinates with respect to S and obtain the manipulated
+ // quadrature points we need. Obviously, the corresponding quadrature points
+ // in world coordinates coincide. Thus the numerical integration on S with the
+ // manipulated quadrature formula gives us the same result as the numerical
+ // integration on S' with the original quadrature formula up to the determinant.
+ // This method for integration on subelements allows the reuse of the routines for
+ // the integration on elements.
+ //
+ // Main routines:
+ // ScalableQuadrature() - Create a new quadrature formula which is a copy of the
+ // original quadrature formula and store the original
+ // quadrature points in \ref oldLambda.
+ // scaleQuadrature() - Manipulate the original quadrature formula for a
+ // subelement.
+ // ===============================================================================
+ class ScalableQuadrature : public Quadrature
+ {
+ public:
+ MEMORY_MANAGED(ScalableQuadrature);
- /**
- * Constructor
- */
- ScalableQuadrature(Quadrature *quadrature);
+ /**
+ * Constructor
+ */
+ ScalableQuadrature(Quadrature *quadrature);
- /**
- * Destructor
- */
- ~ScalableQuadrature()
- {
- DELETE oldLambda;
- }
+ /**
+ * Destructor
+ */
+ ~ScalableQuadrature() {
+ DELETE oldLambda;
+ }
- /**
- * Manipulates the quadrature points for the assemblage of a subelement.
- */
- void scaleQuadrature(const SubElInfo& subElInfo);
+ /**
+ * Manipulates the quadrature points for the assemblage of a subelement.
+ */
+ void scaleQuadrature(const SubElInfo& subElInfo);
+
+ /**
+ * Get b-th coordinate of the a-th original quadrature point.
+ */
+ inline const double getOldLambda(int a, int b) const {
+ if (oldLambda)
+ return (*oldLambda)[a][b];
+ else
+ return 0.0;
+ };
- /**
- * Get b-th coordinate of the a-th original quadrature point.
- */
- inline const double getOldLambda(int a,int b) const {
- if (oldLambda) return (*oldLambda)[a][b];
- else return 0.0;
- };
+ protected:
+ /**
+ * Original quadrature points.
+ */
+ VectorOfFixVecs<DimVec<double> > *oldLambda;
+ };
-protected:
- /**
- * Original quadrature points.
- */
- VectorOfFixVecs<DimVec<double> > *oldLambda;
-};
+}
#endif // AMDIS_SCALABLEQUADRATURE_H
diff --git a/AMDiS/compositeFEM/src/SubElInfo.cc b/AMDiS/compositeFEM/src/SubElInfo.cc
index 876ed8a4fcf80c2a59e7387d9d993f1f780bed54..d0056c3fe3df8607cdd97d93702d7e83cc1200a9 100644
--- a/AMDiS/compositeFEM/src/SubElInfo.cc
+++ b/AMDiS/compositeFEM/src/SubElInfo.cc
@@ -2,34 +2,37 @@
#include "ElInfo.h"
#include "Mesh.h"
-SubElInfo::SubElInfo(VectorOfFixVecs<DimVec<double> > *lambda_,
- const ElInfo *elInfo_)
- : elInfo(elInfo_)
-{
- FUNCNAME("SubElInfo::SubElInfo");
+namespace AMDiS {
- int numPoints = lambda_->getSize();
- int dim = elInfo_->getMesh()->getDim();
- int i;
+ SubElInfo::SubElInfo(VectorOfFixVecs<DimVec<double> > *lambda_,
+ const ElInfo *elInfo_)
+ : elInfo(elInfo_)
+ {
+ FUNCNAME("SubElInfo::SubElInfo");
- TEST_EXIT(numPoints == dim+1)("invalid number of vertices of subelement\n");
+ int nPoints = lambda_->getSize();
+ int dim = elInfo_->getMesh()->getDim();
- FixVec<WorldVector<double>, VERTEX> worldCoords(dim, NO_INIT);
+ TEST_EXIT(nPoints == dim + 1)
+ ("invalid number of vertices of subelement\n");
- lambda = NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
- for (i = 0; i <= dim; i++) {
- (*lambda)[i] = (*lambda_)[i];
- }
+ FixVec<WorldVector<double>, VERTEX> worldCoords(dim, NO_INIT);
- /**
- * Get worldcoordinates of the vertices of subelement in order to
- * calculate the corresponding determinant.
- */
- for (i = 0; i < numPoints; i++) {
- elInfo->coordToWorld((const DimVec<double>)((*lambda)[i]),
- &(worldCoords[i]));
- }
+ lambda = NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
+ for (int i = 0; i <= dim; i++) {
+ (*lambda)[i] = (*lambda_)[i];
+ }
- det = ElInfo::calcDet(worldCoords);
-}
+ /**
+ * Get worldcoordinates of the vertices of subelement in order to
+ * calculate the corresponding determinant.
+ */
+ for (int i = 0; i < nPoints; i++) {
+ elInfo->coordToWorld((const DimVec<double>)((*lambda)[i]),
+ &(worldCoords[i]));
+ }
+ det = ElInfo::calcDet(worldCoords);
+ }
+
+}
diff --git a/AMDiS/compositeFEM/src/SubElInfo.h b/AMDiS/compositeFEM/src/SubElInfo.h
index acacdf360b1536a5a00e40c774b450eebc588d67..dd9706e49dc897e095b84256b6fe481b59be015d 100644
--- a/AMDiS/compositeFEM/src/SubElInfo.h
+++ b/AMDiS/compositeFEM/src/SubElInfo.h
@@ -7,92 +7,91 @@
#include "FixVec.h"
namespace AMDiS {
- class ElInfo;
-}
+ class ElInfo;
-using namespace AMDiS;
-using namespace std;
+ // ============================================================================
+ // ===== class SubElInfo ======================================================
+ // ============================================================================
+ //
+ // Class description:
+ // The class \ref SubElInfo holds all information on a subelement (element to
+ // which it belongs, its vertices in barycentric coordinates with respect to
+ // the element, corresponding determinant, ...).
+ //
+ // Main routines:
+ // SubElInfo() - Creates a \ref SubElInfo for a subelement of the element
+ // \ref elInfo. The vertices of the subelement given in
+ // barycentric coordinates with respect to element are handed over
+ // to this routine. The determinant corresponding to subelement
+ // is calculated.
+ // ============================================================================
+ class SubElInfo
+ {
+ public:
+ MEMORY_MANAGED(SubElInfo);
-// ============================================================================
-// ===== class SubElInfo ======================================================
-// ============================================================================
-//
-// Class description:
-// The class \ref SubElInfo holds all information on a subelement (element to
-// which it belongs, its vertices in barycentric coordinates with respect to
-// the element, corresponding determinant, ...).
-//
-// Main routines:
-// SubElInfo() - Creates a \ref SubElInfo for a subelement of the element
-// \ref elInfo. The vertices of the subelement given in
-// barycentric coordinates with respect to element are handed over
-// to this routine. The determinant corresponding to subelement
-// is calculated.
-// ============================================================================
-class SubElInfo
-{
-public:
- MEMORY_MANAGED(SubElInfo);
+ /**
+ * Constructor
+ */
+ SubElInfo(VectorOfFixVecs<DimVec<double> > *lambda_,
+ const ElInfo *elInfo_);
- /**
- * Constructor
- */
- SubElInfo(VectorOfFixVecs<DimVec<double> > *lambda_,
- const ElInfo *elInfo_);
+ /**
+ * Destructor.
+ */
+ ~SubElInfo() {
+ DELETE lambda;
+ };
- /**
- * Destructor.
- */
- ~SubElInfo()
- {
- DELETE lambda;
- };
- /**
- * Get b-th coordinate of the a-th vertex of subelement (barycentric
- * coordinates).
- */
- inline const double getLambda(int a,int b) const {
- if (lambda) return (*lambda)[a][b];
- else return 0.0;
- };
+ /**
+ * Get b-th coordinate of the a-th vertex of subelement (barycentric
+ * coordinates).
+ */
+ inline const double getLambda(int a,int b) const {
+ if (lambda)
+ return (*lambda)[a][b];
+ else
+ return 0.0;
+ };
- /**
- * Get coordinates of a-th vertex of subelement (barycentric coordinates).
- */
- inline const DimVec<double>& getLambda(int a) const {
- return (*lambda)[a];
- };
+ /**
+ * Get coordinates of a-th vertex of subelement (barycentric coordinates).
+ */
+ inline const DimVec<double>& getLambda(int a) const {
+ return (*lambda)[a];
+ };
- /**
- * Get coordinates of all vertices of subelement (barycentric coordinates).
- */
- inline VectorOfFixVecs<DimVec<double> > *getLambda() const {
- return lambda;
- };
+ /**
+ * Get coordinates of all vertices of subelement (barycentric coordinates).
+ */
+ inline VectorOfFixVecs<DimVec<double> > *getLambda() const {
+ return lambda;
+ };
- /**
- * Get determinant corresponding to subelement.
- */
- inline const double getDet() const {
- return det;
- };
+ /**
+ * Get determinant corresponding to subelement.
+ */
+ inline const double getDet() const {
+ return det;
+ };
-protected:
+ protected:
- /**
- * Contains elInfo of the element that contains subelement.
- */
- const ElInfo *elInfo;
+ /**
+ * Contains elInfo of the element that contains subelement.
+ */
+ const ElInfo *elInfo;
- /**
- * Barycentrc coordinates of the vertices of subelement.
- */
- VectorOfFixVecs<DimVec<double> > *lambda;
+ /**
+ * Barycentrc coordinates of the vertices of subelement.
+ */
+ VectorOfFixVecs<DimVec<double> > *lambda;
- /**
- * Determinant corresponding to the (World-)subelement.
- */
- double det;
-};
+ /**
+ * Determinant corresponding to the (World-)subelement.
+ */
+ double det;
+ };
+}
#endif // AMDIS_SUBELINFO_H
diff --git a/AMDiS/compositeFEM/src/SubElementAssembler.cc b/AMDiS/compositeFEM/src/SubElementAssembler.cc
index 810d45230146f2987d16fcc1bfc8ed143d75214a..883f8c47885e96dbd309ff200e7fd0963c2aed3a 100644
--- a/AMDiS/compositeFEM/src/SubElementAssembler.cc
+++ b/AMDiS/compositeFEM/src/SubElementAssembler.cc
@@ -5,195 +5,184 @@
#include "ScalableQuadrature.h"
#include "SubPolytope.h"
-SubElementAssembler::SubElementAssembler(Operator *op,
- const FiniteElemSpace *rowFESpace_,
- const FiniteElemSpace *colFESpace_)
- : StandardAssembler(op, NULL, NULL, NULL, NULL, rowFESpace_, colFESpace_)
-{
- /**
- * Create a scalable quadrature for subassembler and replace the original
- * quadrature of the subassembler with the scalable quadrature.
- */
-
- if (zeroOrderAssembler) {
- checkQuadratures();
- zeroOrderScalableQuadrature =
- NEW ScalableQuadrature(zeroOrderAssembler->getQuadrature());
- zeroOrderAssembler->setQuadrature(zeroOrderScalableQuadrature);
- }
- else {
- zeroOrderScalableQuadrature = NULL;
- }
+namespace AMDiS {
+
+ SubElementAssembler::SubElementAssembler(Operator *op,
+ const FiniteElemSpace *rowFESpace_,
+ const FiniteElemSpace *colFESpace_)
+ : StandardAssembler(op, NULL, NULL, NULL, NULL, rowFESpace_, colFESpace_)
+ {
+ /**
+ * Create a scalable quadrature for subassembler and replace the original
+ * quadrature of the subassembler with the scalable quadrature.
+ */
- if (firstOrderAssemblerGrdPsi) {
- checkQuadratures();
- firstOrderGrdPsiScalableQuadrature =
- NEW ScalableQuadrature(firstOrderAssemblerGrdPsi->getQuadrature());
- firstOrderAssemblerGrdPsi->setQuadrature(firstOrderGrdPsiScalableQuadrature);
- }
- else {
- firstOrderGrdPsiScalableQuadrature = NULL;
- }
+ if (zeroOrderAssembler) {
+ checkQuadratures();
+ zeroOrderScalableQuadrature =
+ NEW ScalableQuadrature(zeroOrderAssembler->getQuadrature());
+ zeroOrderAssembler->setQuadrature(zeroOrderScalableQuadrature);
+ } else {
+ zeroOrderScalableQuadrature = NULL;
+ }
- if (firstOrderAssemblerGrdPhi) {
- checkQuadratures();
- firstOrderGrdPhiScalableQuadrature =
- NEW ScalableQuadrature(firstOrderAssemblerGrdPhi->getQuadrature());
- firstOrderAssemblerGrdPhi->setQuadrature(firstOrderGrdPhiScalableQuadrature);
- }
- else {
- firstOrderGrdPhiScalableQuadrature = NULL;
- }
+ if (firstOrderAssemblerGrdPsi) {
+ checkQuadratures();
+ firstOrderGrdPsiScalableQuadrature =
+ NEW ScalableQuadrature(firstOrderAssemblerGrdPsi->getQuadrature());
+ firstOrderAssemblerGrdPsi->setQuadrature(firstOrderGrdPsiScalableQuadrature);
+ } else {
+ firstOrderGrdPsiScalableQuadrature = NULL;
+ }
- if (secondOrderAssembler) {
- checkQuadratures();
- secondOrderScalableQuadrature =
- NEW ScalableQuadrature(secondOrderAssembler->getQuadrature());
- secondOrderAssembler->setQuadrature(secondOrderScalableQuadrature);
- }
- else {
- secondOrderScalableQuadrature = NULL;
- }
-}
+ if (firstOrderAssemblerGrdPhi) {
+ checkQuadratures();
+ firstOrderGrdPhiScalableQuadrature =
+ NEW ScalableQuadrature(firstOrderAssemblerGrdPhi->getQuadrature());
+ firstOrderAssemblerGrdPhi->setQuadrature(firstOrderGrdPhiScalableQuadrature);
+ } else {
+ firstOrderGrdPhiScalableQuadrature = NULL;
+ }
-void
-SubElementAssembler::scaleQuadratures(const SubElInfo& subElInfo)
-{
- if (zeroOrderAssembler) {
- zeroOrderScalableQuadrature->scaleQuadrature(subElInfo);
- }
- if (firstOrderAssemblerGrdPsi) {
- firstOrderGrdPsiScalableQuadrature->scaleQuadrature(subElInfo);
- }
- if (firstOrderAssemblerGrdPhi) {
- firstOrderGrdPhiScalableQuadrature->scaleQuadrature(subElInfo);
+ if (secondOrderAssembler) {
+ checkQuadratures();
+ secondOrderScalableQuadrature =
+ NEW ScalableQuadrature(secondOrderAssembler->getQuadrature());
+ secondOrderAssembler->setQuadrature(secondOrderScalableQuadrature);
+ } else {
+ secondOrderScalableQuadrature = NULL;
+ }
}
- if (secondOrderAssembler) {
- secondOrderScalableQuadrature->scaleQuadrature(subElInfo);
+
+ void SubElementAssembler::scaleQuadratures(const SubElInfo& subElInfo)
+ {
+ if (zeroOrderAssembler) {
+ zeroOrderScalableQuadrature->scaleQuadrature(subElInfo);
+ }
+ if (firstOrderAssemblerGrdPsi) {
+ firstOrderGrdPsiScalableQuadrature->scaleQuadrature(subElInfo);
+ }
+ if (firstOrderAssemblerGrdPhi) {
+ firstOrderGrdPhiScalableQuadrature->scaleQuadrature(subElInfo);
+ }
+ if (secondOrderAssembler) {
+ secondOrderScalableQuadrature->scaleQuadrature(subElInfo);
+ }
}
-}
-void
-SubElementAssembler::getSubElementVector(SubElInfo *subElInfo,
- const ElInfo *elInfo,
- ElementVector *userVec)
-{
- double corrFactor;
-
- /**
- * Manipulate the quadratures of the SubAssemblers for subelement.
- */
- scaleQuadratures(*subElInfo);
-
- calculateElementVector(elInfo, userVec);
-
- /**
- * The integration has been performed with a quadrature living on element. The
- * determinant of element has been used instead of the determinant of subelement. Thus
- * the result must be corrected with respect to subelement.
- */
- corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
- int i;
- for ( i = 0; i < nRow; i++) {
- (*userVec)[i] *= corrFactor;
+ void SubElementAssembler::getSubElementVector(SubElInfo *subElInfo,
+ const ElInfo *elInfo,
+ ElementVector *userVec)
+ {
+ double corrFactor;
+
+ /**
+ * Manipulate the quadratures of the SubAssemblers for subelement.
+ */
+ scaleQuadratures(*subElInfo);
+
+ calculateElementVector(elInfo, userVec);
+
+ /**
+ * The integration has been performed with a quadrature living on element. The
+ * determinant of element has been used instead of the determinant of subelement. Thus
+ * the result must be corrected with respect to subelement.
+ */
+ corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
+ for (int i = 0; i < nRow; i++) {
+ (*userVec)[i] *= corrFactor;
+ }
}
-}
-void
-SubElementAssembler::getSubElementMatrix(SubElInfo *subElInfo,
- const ElInfo *elInfo,
- ElementMatrix *userMat)
-{
- double corrFactor;
-
- /**
- * Manipulate the quadratures of the SubAssemblers for subelement.
- */
- scaleQuadratures(*subElInfo);
-
- /**
- * Integrate using the manipulated quadrature.
- */
- calculateElementMatrix(elInfo, userMat);
-
- /**
- * The integration has been performed with a quadrature living on element. The
- * determinant of element has been used instead of the determinant of subelement.
- * Thus the result must be corrected with respect to subelement.
- */
- corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
- int i, j;
- for ( i = 0; i < nRow; i++) {
- for ( j = 0; j < nCol; j++) {
- (*userMat)[i][j] *= corrFactor;
+ void SubElementAssembler::getSubElementMatrix(SubElInfo *subElInfo,
+ const ElInfo *elInfo,
+ ElementMatrix *userMat)
+ {
+ /**
+ * Manipulate the quadratures of the SubAssemblers for subelement.
+ */
+ scaleQuadratures(*subElInfo);
+
+ /**
+ * Integrate using the manipulated quadrature.
+ */
+ calculateElementMatrix(elInfo, userMat);
+
+ /**
+ * The integration has been performed with a quadrature living on element. The
+ * determinant of element has been used instead of the determinant of subelement.
+ * Thus the result must be corrected with respect to subelement.
+ */
+ double corrFactor = subElInfo->getDet() / fabs(elInfo->getDet());
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++) {
+ (*userMat)[i][j] *= corrFactor;
+ }
}
}
-}
-void
-SubElementAssembler::getSubPolytopeVector(SubPolytope *subPolytope,
- SubElementAssembler *subElementAssembler,
- const ElInfo *elInfo,
- ElementVector *subPolVec)
-{
- /**
- * Note: There is no reset of subPolVec.
- */
- std::vector<SubElInfo *>::iterator it;
- ElementVector *subElVec = NEW ElementVector(nRow);
- int i;
-
- /**
- * Assemble for each subelement of subpolytope.
- */
- for (it = subPolytope->getSubElementsBegin();
- it != subPolytope->getSubElementsEnd();
- it++) {
- subElVec->set(0.0);
- subElementAssembler->getSubElementVector(*it, elInfo, subElVec);
+ void SubElementAssembler::getSubPolytopeVector(SubPolytope *subPolytope,
+ SubElementAssembler *subElementAssembler,
+ const ElInfo *elInfo,
+ ElementVector *subPolVec)
+ {
+ /**
+ * Note: There is no reset of subPolVec.
+ */
+ std::vector<SubElInfo *>::iterator it;
+ ElementVector *subElVec = NEW ElementVector(nRow);
/**
- * Add results for subelement to total result for subpolytope.
+ * Assemble for each subelement of subpolytope.
*/
- for(i = 0; i < nRow; i++) {
- (*subPolVec)[i] += (*subElVec)[i];
+ for (it = subPolytope->getSubElementsBegin();
+ it != subPolytope->getSubElementsEnd();
+ it++) {
+ subElVec->set(0.0);
+ subElementAssembler->getSubElementVector(*it, elInfo, subElVec);
+
+ /**
+ * Add results for subelement to total result for subpolytope.
+ */
+ for (int i = 0; i < nRow; i++) {
+ (*subPolVec)[i] += (*subElVec)[i];
+ }
}
- }
- DELETE subElVec;
-}
+ DELETE subElVec;
+ }
-void
-SubElementAssembler::getSubPolytopeMatrix(SubPolytope *subPolytope,
- SubElementAssembler *subElementAssembler,
- const ElInfo *elInfo,
- ElementMatrix *subPolMat)
-{
- /**
- * Note: There is no reset of subPolMat.
- */
- std::vector<SubElInfo *>::iterator it;
- ElementMatrix *subElMat = NEW ElementMatrix(nRow, nCol);
- int i,j;
-
- /**
- * Assemble for each subelement of subpolytope.
- */
- for (it = subPolytope->getSubElementsBegin();
- it != subPolytope->getSubElementsEnd();
- it++) {
- subElMat->set(0.0);
- subElementAssembler->getSubElementMatrix(*it, elInfo, subElMat);
+ void SubElementAssembler::getSubPolytopeMatrix(SubPolytope *subPolytope,
+ SubElementAssembler *subElementAssembler,
+ const ElInfo *elInfo,
+ ElementMatrix *subPolMat)
+ {
+ /**
+ * Note: There is no reset of subPolMat.
+ */
+ std::vector<SubElInfo *>::iterator it;
+ ElementMatrix *subElMat = NEW ElementMatrix(nRow, nCol);
/**
- * Add results for subelement to total result for subpolytope.
+ * Assemble for each subelement of subpolytope.
*/
- for(i = 0; i < nRow; i++) {
- for(j = 0; j < nCol; j++) {
- (*subPolMat)[i][j] += (*subElMat)[i][j];
+ for (it = subPolytope->getSubElementsBegin();
+ it != subPolytope->getSubElementsEnd();
+ it++) {
+ subElMat->set(0.0);
+ subElementAssembler->getSubElementMatrix(*it, elInfo, subElMat);
+
+ /**
+ * Add results for subelement to total result for subpolytope.
+ */
+ for (int i = 0; i < nRow; i++) {
+ for (int j = 0; j < nCol; j++) {
+ (*subPolMat)[i][j] += (*subElMat)[i][j];
+ }
}
}
+
+ DELETE subElMat;
}
- DELETE subElMat;
}
diff --git a/AMDiS/compositeFEM/src/SubElementAssembler.h b/AMDiS/compositeFEM/src/SubElementAssembler.h
index 2d995e127fb171d42da7ede840ad1fab489dd5f1..e5697592d80a36395ff2497f5dea00c0ef113bbb 100644
--- a/AMDiS/compositeFEM/src/SubElementAssembler.h
+++ b/AMDiS/compositeFEM/src/SubElementAssembler.h
@@ -6,99 +6,99 @@
#include "MemoryManager.h"
#include "Assembler.h"
#include "SubElInfo.h"
-
#include "ScalableQuadrature.h"
-class SubPolytope;
+namespace AMDiS {
-using namespace AMDiS;
-using namespace std;
+ class SubPolytope;
-// ============================================================================
-// ===== class SubElementAssembler ============================================
-// ============================================================================
-//
-// Class Desription:
-// The class \ref SubElementAssembler holds the routines for the assemblage on
-// subpolytopes and subelements.
-// The integration on a subpolytope takes place by integrating on the
-// subelements and afterwards summing up the results.
-// If S' is a sublement and S the element containing S', the numerical integration
-// on S' is done by integration on the element S with a manipulated quadrature
-// formula and multiplication of the result with a correction term consisting of
-// the determinants corresponding to S and S'. That means, the quadrature points
-// are manipulated in the following way:
-// The original quadrature formula holds quadrature points which are given in
-// barycentric coordinates with respect to S'. Now we express these quadrature
-// points in barycentric coordinates with respect to S and obtain the manipulated
-// quadrature points we need. Obviously, the corresponding quadrature points
-// in world coordinates coincide. Thus the numerical integration on S with the
-// manipulated quadrature formula gives us the same result as the numerical
-// integration on S' with the original quadrature formula up to the determinant.
-// This method for integration on subelements allows the reuse of the routines for
-// the integration on elements.
-//
-// The manipulation of the quadrature formula takes place for the corresponding
-// assembler type (ZeroOrderAssembler, FirstOrderAssemblerGrdPsi, ...).
-//
-// Main routines:
-// SubElementAssembler() - Creates a scalable quadrature for the appropriate
-// assembler type and assigns this quadrature to the
-// assembler.
-// scaleQuadratures() - Manipulates the scalable quadrature of the
-// appropriate assembler type with respect to a
-// subelement.
-// getSubPolytopeVector() - Calculates the righthandside vector for a polytope.
-// getSubPolytopeMatrix() - Calculates the system matrix for a polytope.
-// getSubElementVector() - Calculates the righthandside vector for a subelement.
-// getSubElementMatrix() - Calculates the system matrix for a subelement.
-// ============================================================================
-class SubElementAssembler : public StandardAssembler
-{
- public:
- MEMORY_MANAGED(SubElementAssembler);
+ // ============================================================================
+ // ===== class SubElementAssembler ============================================
+ // ============================================================================
+ //
+ // Class Desription:
+ // The class \ref SubElementAssembler holds the routines for the assemblage on
+ // subpolytopes and subelements.
+ // The integration on a subpolytope takes place by integrating on the
+ // subelements and afterwards summing up the results.
+ // If S' is a sublement and S the element containing S', the numerical integration
+ // on S' is done by integration on the element S with a manipulated quadrature
+ // formula and multiplication of the result with a correction term consisting of
+ // the determinants corresponding to S and S'. That means, the quadrature points
+ // are manipulated in the following way:
+ // The original quadrature formula holds quadrature points which are given in
+ // barycentric coordinates with respect to S'. Now we express these quadrature
+ // points in barycentric coordinates with respect to S and obtain the manipulated
+ // quadrature points we need. Obviously, the corresponding quadrature points
+ // in world coordinates coincide. Thus the numerical integration on S with the
+ // manipulated quadrature formula gives us the same result as the numerical
+ // integration on S' with the original quadrature formula up to the determinant.
+ // This method for integration on subelements allows the reuse of the routines for
+ // the integration on elements.
+ //
+ // The manipulation of the quadrature formula takes place for the corresponding
+ // assembler type (ZeroOrderAssembler, FirstOrderAssemblerGrdPsi, ...).
+ //
+ // Main routines:
+ // SubElementAssembler() - Creates a scalable quadrature for the appropriate
+ // assembler type and assigns this quadrature to the
+ // assembler.
+ // scaleQuadratures() - Manipulates the scalable quadrature of the
+ // appropriate assembler type with respect to a
+ // subelement.
+ // getSubPolytopeVector() - Calculates the righthandside vector for a polytope.
+ // getSubPolytopeMatrix() - Calculates the system matrix for a polytope.
+ // getSubElementVector() - Calculates the righthandside vector for a subelement.
+ // getSubElementMatrix() - Calculates the system matrix for a subelement.
+ // ============================================================================
+ class SubElementAssembler : public StandardAssembler
+ {
+ public:
+ MEMORY_MANAGED(SubElementAssembler);
- SubElementAssembler(Operator *op,
- const FiniteElemSpace *rowFESpace_,
- const FiniteElemSpace *colFESpace_=NULL);
+ SubElementAssembler(Operator *op,
+ const FiniteElemSpace *rowFESpace_,
+ const FiniteElemSpace *colFESpace_=NULL);
- virtual ~SubElementAssembler()
- {
- if (zeroOrderScalableQuadrature)
- DELETE zeroOrderScalableQuadrature;
- if (firstOrderGrdPsiScalableQuadrature)
- DELETE firstOrderGrdPsiScalableQuadrature;
- if (firstOrderGrdPhiScalableQuadrature)
- DELETE firstOrderGrdPhiScalableQuadrature;
- if (secondOrderScalableQuadrature)
- DELETE secondOrderScalableQuadrature;
- };
+ virtual ~SubElementAssembler()
+ {
+ if (zeroOrderScalableQuadrature)
+ DELETE zeroOrderScalableQuadrature;
+ if (firstOrderGrdPsiScalableQuadrature)
+ DELETE firstOrderGrdPsiScalableQuadrature;
+ if (firstOrderGrdPhiScalableQuadrature)
+ DELETE firstOrderGrdPhiScalableQuadrature;
+ if (secondOrderScalableQuadrature)
+ DELETE secondOrderScalableQuadrature;
+ };
+
+ void scaleQuadratures(const SubElInfo& subElInfo);
- void scaleQuadratures(const SubElInfo& subElInfo);
+ void getSubElementVector(SubElInfo *subElInfo,
+ const ElInfo *elInfo,
+ ElementVector *userVec);
- void getSubElementVector(SubElInfo *subElInfo,
- const ElInfo *elInfo,
- ElementVector *userVec);
+ void getSubElementMatrix(SubElInfo *subElInfo,
+ const ElInfo *elInfo,
+ ElementMatrix *userMat);
- void getSubElementMatrix(SubElInfo *subElInfo,
- const ElInfo *elInfo,
- ElementMatrix *userMat);
+ void getSubPolytopeVector(SubPolytope *subPolytope,
+ SubElementAssembler *subElementAssembler,
+ const ElInfo *elInfo,
+ ElementVector *userVec);
- void getSubPolytopeVector(SubPolytope *subPolytope,
- SubElementAssembler *subElementAssembler,
- const ElInfo *elInfo,
- ElementVector *userVec);
+ void getSubPolytopeMatrix(SubPolytope *subPolytope,
+ SubElementAssembler *subElementAssembler,
+ const ElInfo *elInfo,
+ ElementMatrix *userMat);
- void getSubPolytopeMatrix(SubPolytope *subPolytope,
- SubElementAssembler *subElementAssembler,
- const ElInfo *elInfo,
- ElementMatrix *userMat);
+ protected:
+ ScalableQuadrature *zeroOrderScalableQuadrature;
+ ScalableQuadrature *firstOrderGrdPsiScalableQuadrature;
+ ScalableQuadrature *firstOrderGrdPhiScalableQuadrature;
+ ScalableQuadrature *secondOrderScalableQuadrature;
+ };
-protected:
- ScalableQuadrature *zeroOrderScalableQuadrature;
- ScalableQuadrature *firstOrderGrdPsiScalableQuadrature;
- ScalableQuadrature *firstOrderGrdPhiScalableQuadrature;
- ScalableQuadrature *secondOrderScalableQuadrature;
-};
+}
#endif // AMDIS_SUBELEMENTASSEMBLER_H
diff --git a/AMDiS/compositeFEM/src/SubPolytope.cc b/AMDiS/compositeFEM/src/SubPolytope.cc
index 5f9c8de7ce22393ad645b0804effc32de3f07606..4403dc7f6319880ea0e86c1e37a3d49e0fa9aaaa 100644
--- a/AMDiS/compositeFEM/src/SubPolytope.cc
+++ b/AMDiS/compositeFEM/src/SubPolytope.cc
@@ -3,74 +3,76 @@
#include "SubElInfo.h"
#include "SubPolytope.h"
-bool SubPolytope::checkIntPoints()
-{
- ////////////////////////////////////////////////////////////////////////////
- //
- // Do the points in intPoints lie on edges ? And are they inner points,
- // i.e. they aren't vertices ?
- //
- // Return value: true - yes
- // false - no
- ////////////////////////////////////////////////////////////////////////////
-
- int zeroCounter;
-
- for (int i = 0; i < numIntPoints; i++) {
- zeroCounter = 0;
- for (int j = 0; j < dim+1; j++) {
- if (fabs((*intPoints)[i][j]) <= 1.e-15 ) {
- zeroCounter++;
+namespace AMDiS {
+
+ bool SubPolytope::checkIntPoints()
+ {
+ ////////////////////////////////////////////////////////////////////////////
+ //
+ // Do the points in intPoints lie on edges ? And are they inner points,
+ // i.e. they aren't vertices ?
+ //
+ // Return value: true - yes
+ // false - no
+ ////////////////////////////////////////////////////////////////////////////
+
+ int zeroCounter;
+
+ for (int i = 0; i < numIntPoints; i++) {
+ zeroCounter = 0;
+ for (int j = 0; j < dim+1; j++) {
+ if (fabs((*intPoints)[i][j]) <= 1.e-15 ) {
+ zeroCounter++;
+ }
}
- }
- /**
- * Dimension 1
- *
- * "Inner" points on edges aren't equal to 0.0 in any component.
- */
- if (dim == 1 && zeroCounter != 0 ) {
- return false;
- }
+ /**
+ * Dimension 1
+ *
+ * "Inner" points on edges aren't equal to 0.0 in any component.
+ */
+ if (dim == 1 && zeroCounter != 0 ) {
+ return false;
+ }
- /**
- * Dimension 2
- *
- * "Inner" points on edges are equal to 0.0 in exactly 1 component.
- */
- if (dim == 2 && zeroCounter != 1) {
- return false;
- }
+ /**
+ * Dimension 2
+ *
+ * "Inner" points on edges are equal to 0.0 in exactly 1 component.
+ */
+ if (dim == 2 && zeroCounter != 1) {
+ return false;
+ }
- /**
- * Dimension 3
- *
- * "Inner" points on edges are equal to 0.0 in exactly 2 components.
- */
- if (dim == 3 && zeroCounter != 2) {
- return false;
+ /**
+ * Dimension 3
+ *
+ * "Inner" points on edges are equal to 0.0 in exactly 2 components.
+ */
+ if (dim == 3 && zeroCounter != 2) {
+ return false;
+ }
}
- }
- return true;
-}
+ return true;
+ }
-SubPolytope::SubPolytope(const ElInfo *elInfo_,
- VectorOfFixVecs<DimVec<double> > *intPoints_,
- int numIntPoints_,
- const int &indexElVertInPol_)
- : elInfo(elInfo_),
- intPoints(intPoints_),
- numIntPoints(numIntPoints_)
-{
- FUNCNAME("SubPolytope::SubPolytope()");
-
- dim = (*intPoints_)[0].getSize() - 1;
-
- TEST_EXIT((dim == 1 && numIntPoints == 1) ||
- (dim == 2 && numIntPoints == 2) ||
- (dim == 3 && numIntPoints == 3) ||
- (dim == 3 && numIntPoints == 4))
+ SubPolytope::SubPolytope(const ElInfo *elInfo_,
+ VectorOfFixVecs<DimVec<double> > *intPoints_,
+ int numIntPoints_,
+ const int &indexElVertInPol_)
+ : elInfo(elInfo_),
+ intPoints(intPoints_),
+ numIntPoints(numIntPoints_)
+ {
+ FUNCNAME("SubPolytope::SubPolytope()");
+
+ dim = (*intPoints_)[0].getSize() - 1;
+
+ TEST_EXIT((dim == 1 && numIntPoints == 1) ||
+ (dim == 2 && numIntPoints == 2) ||
+ (dim == 3 && numIntPoints == 3) ||
+ (dim == 3 && numIntPoints == 4))
("invalid number of intersection points\n");
/**
@@ -101,389 +103,391 @@ SubPolytope::SubPolytope(const ElInfo *elInfo_,
ERROR_EXIT("invalid dimension\n");
break;
}
-}
+ }
-void SubPolytope::createSubElementPolytopeIsSubElement1D(int indexElVertInPol)
-{
- //**************************************************************************
- // The intersection of the one-dimensional element (interval) divides
- // element into two subelements. indexElVertInPol indicates which
- // subelement to take.
- //**************************************************************************
+ void SubPolytope::createSubElementPolytopeIsSubElement1D(int indexElVertInPol)
+ {
+ //**************************************************************************
+ // The intersection of the one-dimensional element (interval) divides
+ // element into two subelements. indexElVertInPol indicates which
+ // subelement to take.
+ //**************************************************************************
- FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement1D()");
+ FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement1D()");
- TEST_EXIT(dim == 1 && numIntPoints == 1)("invalid call of this routine\n");
+ TEST_EXIT(dim == 1 && numIntPoints == 1)("invalid call of this routine\n");
- VectorOfFixVecs<DimVec<double> > *subElVertices =
- NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
- DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
+ VectorOfFixVecs<DimVec<double> > *subElVertices =
+ NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
+ DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
- /**
- * Get the vertex which - with the intersection point in intPoints - forms
- * a subelement of element.
- */
- if (indexElVertInPol == 0) {
/**
- * The vertex in element with barycentric coordinates (1,0) is in
- * subelement.
+ * Get the vertex which - with the intersection point in intPoints - forms
+ * a subelement of element.
*/
- vertex[1] = 0.0;
- }
- else {
+ if (indexElVertInPol == 0) {
+ /**
+ * The vertex in element with barycentric coordinates (1,0) is in
+ * subelement.
+ */
+ vertex[1] = 0.0;
+ }
+ else {
+ /**
+ * The vertex in element with barycentric coordinates (0,1) is in
+ * subelement.
+ */
+ vertex[0] = 0.0;
+ }
+
/**
- * The vertex in element with barycentric coordinates (0,1) is in
- * subelement.
+ * Collect the vertices of the subelement in subElVertices.
+ *
+ * Note: The routines CompositeFEMOperator::getElementMatrix and
+ * CompositeFEMOperator::getElementVector expect the first vertice
+ * of a subelement to be a vertice of the corresponding element and
+ * not to be an intersection point.
*/
- vertex[0] = 0.0;
- }
+ (*subElVertices)[0] = vertex;
+ (*subElVertices)[dim] = (*intPoints)[0];
- /**
- * Collect the vertices of the subelement in subElVertices.
- *
- * Note: The routines CompositeFEMOperator::getElementMatrix and
- * CompositeFEMOperator::getElementVector expect the first vertice
- * of a subelement to be a vertice of the corresponding element and
- * not to be an intersection point.
- */
- (*subElVertices)[0] = vertex;
- (*subElVertices)[dim] = (*intPoints)[0];
+ /**
+ * Create a new ElInfo for the subelement.
+ */
+ subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
- /**
- * Create a new ElInfo for the subelement.
- */
- subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
+ TEST_EXIT( subElements.size() == 1 )("error in creating subelements");
+ numSubElements = 1;
- TEST_EXIT( subElements.size() == 1 )("error in creating subelements");
- numSubElements = 1;
+ DELETE subElVertices;
+ }
- DELETE subElVertices;
-}
+ void SubPolytope::createSubElementPolytopeIsSubElement2D3D()
+ {
+ //**************************************************************************
+ // The intersection with element produced dim intersection points.
+ // Thus there is exactly one vertice in element which - with the
+ // intersection points - forms a subelement of element.
+ // This routine determines this vertex and creates a new object of SubElInfo
+ // for the subelement stored in subElements.
+ //
+ // How does it work ?
+ // All intersection points lie on edges of element and aren't vertices. The
+ // missing vertex is the intersection point of these edges.
+ // Lying on edges means, that each intersection point has exactly one
+ // component equal to 0.0. The missing vertex is equal to 0.0 in all these
+ // components.
+ //**************************************************************************
+
+ FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement2D3D()");
+
+ TEST_EXIT((dim == 2 && numIntPoints == 2) ||
+ (dim == 3 && numIntPoints == 3))
+ ("invalid call of this routine\n");
+
+ VectorOfFixVecs<DimVec<double> >*subElVertices =
+ NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
+ DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
-void SubPolytope::createSubElementPolytopeIsSubElement2D3D()
-{
- //**************************************************************************
- // The intersection with element produced dim intersection points.
- // Thus there is exactly one vertice in element which - with the
- // intersection points - forms a subelement of element.
- // This routine determines this vertex and creates a new object of SubElInfo
- // for the subelement stored in subElements.
- //
- // How does it work ?
- // All intersection points lie on edges of element and aren't vertices. The
- // missing vertex is the intersection point of these edges.
- // Lying on edges means, that each intersection point has exactly one
- // component equal to 0.0. The missing vertex is equal to 0.0 in all these
- // components.
- //**************************************************************************
-
- FUNCNAME("SubPolytope::createSubElementPolytopeIsSubElement2D3D()");
-
- TEST_EXIT((dim == 2 && numIntPoints == 2) ||
- (dim == 3 && numIntPoints == 3))
- ("invalid call of this routine\n");
-
- VectorOfFixVecs<DimVec<double> >*subElVertices =
- NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
- DimVec<double> vertex(dim, DEFAULT_VALUE, 1.0);
-
- /**
- * Get the vertex which - with the intersection points intPoints - forms
- * a subelement of element.
- */
- for (int i = 0; i < numIntPoints; i++) {
- for (int j = 0; j < dim+1; j++) {
- if ( fabs((*intPoints)[i][j]) <= 1.e-15 ) {
- vertex[j] = 0.0;
- };
+ /**
+ * Get the vertex which - with the intersection points intPoints - forms
+ * a subelement of element.
+ */
+ for (int i = 0; i < numIntPoints; i++) {
+ for (int j = 0; j < dim+1; j++) {
+ if ( fabs((*intPoints)[i][j]) <= 1.e-15 ) {
+ vertex[j] = 0.0;
+ };
+ }
}
- }
- /**
- * Collect the vertices of the subelement in subElVertices.
- *
- * Note: The routines CompositeFEMOperator::getElementMatrix and
- * CompositeFEMOperator::getElementVector expect the first vertice
- * of a subelement to be a vertice of the corresponding element and
- * not to be an intersection point.
- */
- (*subElVertices)[0] = vertex;
- for (int i = 0; i < numIntPoints; i++) {
- (*subElVertices)[i+1] = (*intPoints)[i];
- }
+ /**
+ * Collect the vertices of the subelement in subElVertices.
+ *
+ * Note: The routines CompositeFEMOperator::getElementMatrix and
+ * CompositeFEMOperator::getElementVector expect the first vertice
+ * of a subelement to be a vertice of the corresponding element and
+ * not to be an intersection point.
+ */
+ (*subElVertices)[0] = vertex;
+ for (int i = 0; i < numIntPoints; i++) {
+ (*subElVertices)[i+1] = (*intPoints)[i];
+ }
- /**
- * Create a new ElInfo for the subelement.
- */
- subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
+ /**
+ * Create a new ElInfo for the subelement.
+ */
+ subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
- TEST_EXIT( subElements.size() == 1 )("error in creating subelements");
- numSubElements = 1;
+ TEST_EXIT( subElements.size() == 1 )("error in creating subelements");
+ numSubElements = 1;
- DELETE subElVertices;
-}
+ DELETE subElVertices;
+ }
-int SubPolytope::getIndexSecondFaceIntPoint0(int indexFirstFace, int dim)
-{
- for (int i = 0; i < dim+1; i++) {
- if ( fabs((*intPoints)[0][i]) <= 1.e-15 && i != indexFirstFace ) {
- return i;
+ int SubPolytope::getIndexSecondFaceIntPoint0(int indexFirstFace, int dim)
+ {
+ for (int i = 0; i < dim+1; i++) {
+ if ( fabs((*intPoints)[0][i]) <= 1.e-15 && i != indexFirstFace ) {
+ return i;
+ }
}
- }
- ERROR_EXIT("couldn't determine the second face for IntPoint0\n");
- return -1;
-}
+ ERROR_EXIT("couldn't determine the second face for IntPoint0\n");
+ return -1;
+ }
-void SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
-{
- //**************************************************************************
- // The intersection with element produced four intersection points. Thus the
- // intersection doesn't induce a subelement. This routine divides the
- // subpolytope given by the intersection into three subelements.
- //
- // How does it work ?
- // The intersection points and two vertices of element form a subpolytope
- // of element. First of all, we determine these vertices, and call them
- // A and B. Then we sort the intersection points (S_0, S_1, S_2, S_3)
- // in the following way:
- // S_0 is the first intersection point in the creation-list /ref intPoints_.
- // A and S_0 lie in the face of element opposite to B.
- // S_0 and S_1 are neighbours of A.
- // S_2 is opposite to S_0 in the intersection plane.
- // S_1 and S_3 are neighbours of A.
- // Then the subelements are:
- // A - B - S_0 - S_1 , B - S_0 - S_1 - S_2 , B - S_0 - S_2 - S_3 .
- //
- // The index of one vertex of element that is in subpolytope is handed to
- // this routine. The subpolytope is determined uniquely by this vertex and
- // the intersection points.
- //**************************************************************************
-
- FUNCNAME("SubPolytope::createSubElementForPolytope3D");
-
- TEST_EXIT(dim == 3 && numIntPoints == 4)("invalid call of this routine\n");
-
- TEST_EXIT(0 <= indexElVertInPol1 && indexElVertInPol1 <= 3)
- ("invalid index for vertex of a tetrahedron");
-
- VectorOfFixVecs<DimVec<double> > *subElVertices =
- NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
- DimVec<double> vertexA(dim, DEFAULT_VALUE, 0.0);
- DimVec<double> vertexB(dim, DEFAULT_VALUE, 0.0);
-
- int indexElVertInPol2; // index of second vertex of element lying in
- // subpolytope
- // The vertices of element (3D -> tetrahedron) are
- // indexed as usual:
- // 0: (1,0,0,0), 1: (0,1,0,0), 2: (0,0,1,0),
- // 3: (0,0,0,1)
- bool intPointOnEdge[4][4] = {{false, false, false, false},
- {false, false, false, false},
- {false, false, false, false},
- {false, false, false, false}};
- // /ref intPointOnEdge[i][j] indicates whether there
- // is an intersection point on the edge from
- // vertice i to vertice j :
- // false : no intersection point
- // true: there is an intersection point
- int indexEdge[2]; // For a vertex lying on an edge indexEdge
- // stores the indices of the two barycentric
- // coordinates which are not equal to zero.
- int indexA = 0;
- int indexB = 0;
- int indexS_0 = 0;
- int indexS_1 = 0;
- int indexS_2 = 0;
- int indexS_3 = 0;
- int indexSecondFaceIntPoint0 = 0;
-
- /**
- * Get the second vertex of element lying in subpolytope.
- *
- * There is exactly one vertex of element which - with vertex
- * indexElVertInPol1 and the intersection points intPoints -
- * forms a subpolytope of element. It is the vertex adjacent with
- * indexElVertInPol1 whose common edge with indexElVertInPol1
- * doesn't contain an intersection point.
- */
-
- // Get the edges including the intersection points.
- for (int i = 0; i < numIntPoints; i++) {
- int k = 0;
- for (int j = 0; j < dim+1; j++) {
- if (fabs((*intPoints)[i][j]) > 1.e-15 ) {
- indexEdge[k] = j;
- k++;
+ void SubPolytope::createSubElementsForPolytope3D(int indexElVertInPol1)
+ {
+ //**************************************************************************
+ // The intersection with element produced four intersection points. Thus the
+ // intersection doesn't induce a subelement. This routine divides the
+ // subpolytope given by the intersection into three subelements.
+ //
+ // How does it work ?
+ // The intersection points and two vertices of element form a subpolytope
+ // of element. First of all, we determine these vertices, and call them
+ // A and B. Then we sort the intersection points (S_0, S_1, S_2, S_3)
+ // in the following way:
+ // S_0 is the first intersection point in the creation-list /ref intPoints_.
+ // A and S_0 lie in the face of element opposite to B.
+ // S_0 and S_1 are neighbours of A.
+ // S_2 is opposite to S_0 in the intersection plane.
+ // S_1 and S_3 are neighbours of A.
+ // Then the subelements are:
+ // A - B - S_0 - S_1 , B - S_0 - S_1 - S_2 , B - S_0 - S_2 - S_3 .
+ //
+ // The index of one vertex of element that is in subpolytope is handed to
+ // this routine. The subpolytope is determined uniquely by this vertex and
+ // the intersection points.
+ //**************************************************************************
+
+ FUNCNAME("SubPolytope::createSubElementForPolytope3D");
+
+ TEST_EXIT(dim == 3 && numIntPoints == 4)("invalid call of this routine\n");
+
+ TEST_EXIT(0 <= indexElVertInPol1 && indexElVertInPol1 <= 3)
+ ("invalid index for vertex of a tetrahedron");
+
+ VectorOfFixVecs<DimVec<double> > *subElVertices =
+ NEW VectorOfFixVecs<DimVec<double> >(dim, dim+1, NO_INIT);
+ DimVec<double> vertexA(dim, DEFAULT_VALUE, 0.0);
+ DimVec<double> vertexB(dim, DEFAULT_VALUE, 0.0);
+
+ int indexElVertInPol2 = 0; // index of second vertex of element lying in
+ // subpolytope
+ // The vertices of element (3D -> tetrahedron) are
+ // indexed as usual:
+ // 0: (1,0,0,0), 1: (0,1,0,0), 2: (0,0,1,0),
+ // 3: (0,0,0,1)
+ bool intPointOnEdge[4][4] = {{false, false, false, false},
+ {false, false, false, false},
+ {false, false, false, false},
+ {false, false, false, false}};
+ // /ref intPointOnEdge[i][j] indicates whether there
+ // is an intersection point on the edge from
+ // vertice i to vertice j :
+ // false : no intersection point
+ // true: there is an intersection point
+ int indexEdge[2]; // For a vertex lying on an edge indexEdge
+ // stores the indices of the two barycentric
+ // coordinates which are not equal to zero.
+ int indexA = 0;
+ int indexB = 0;
+ int indexS_0 = 0;
+ int indexS_1 = 0;
+ int indexS_2 = 0;
+ int indexS_3 = 0;
+ int indexSecondFaceIntPoint0 = 0;
+
+ /**
+ * Get the second vertex of element lying in subpolytope.
+ *
+ * There is exactly one vertex of element which - with vertex
+ * indexElVertInPol1 and the intersection points intPoints -
+ * forms a subpolytope of element. It is the vertex adjacent with
+ * indexElVertInPol1 whose common edge with indexElVertInPol1
+ * doesn't contain an intersection point.
+ */
+
+ // Get the edges including the intersection points.
+ for (int i = 0; i < numIntPoints; i++) {
+ int k = 0;
+ for (int j = 0; j < dim+1; j++) {
+ if (fabs((*intPoints)[i][j]) > 1.e-15 ) {
+ indexEdge[k] = j;
+ k++;
+ }
}
+ intPointOnEdge[indexEdge[0]][indexEdge[1]] = true;
+ intPointOnEdge[indexEdge[1]][indexEdge[0]] = true;
}
- intPointOnEdge[indexEdge[0]][indexEdge[1]] = true;
- intPointOnEdge[indexEdge[1]][indexEdge[0]] = true;
- }
- // Get the vertex of element adjacent with indexElVertInPol1 whose
- // common edge with indexElVertInPol1 doesn't contain an
- // intersection point, and store it in indexElVertInPol2.
- for (int i = 0; i < dim+1; i++) {
- if (intPointOnEdge[indexElVertInPol1][i] == false &&
- i != indexElVertInPol1 ) {
- indexElVertInPol2 = i;
- break;
+ // Get the vertex of element adjacent with indexElVertInPol1 whose
+ // common edge with indexElVertInPol1 doesn't contain an
+ // intersection point, and store it in indexElVertInPol2.
+ for (int i = 0; i < dim+1; i++) {
+ if (intPointOnEdge[indexElVertInPol1][i] == false &&
+ i != indexElVertInPol1 ) {
+ indexElVertInPol2 = i;
+ break;
+ }
}
- }
- /**
- * Determine A and B, so that intPoint0 is a neighbour of A.
- *
- * In the subpolytope A and two intersection points lie in the face
- * opposite to B. And B and the other two intersection points lie in the
- * face opposite to A.
- */
+ /**
+ * Determine A and B, so that intPoint0 is a neighbour of A.
+ *
+ * In the subpolytope A and two intersection points lie in the face
+ * opposite to B. And B and the other two intersection points lie in the
+ * face opposite to A.
+ */
- if (fabs((*intPoints)[0][indexElVertInPol1]) <= 1.e-15) {
+ if (fabs((*intPoints)[0][indexElVertInPol1]) <= 1.e-15) {
- // (*intPoints)[0] lies in the face opposite to vertex
- // /ref indexElVertInPol1.
+ // (*intPoints)[0] lies in the face opposite to vertex
+ // /ref indexElVertInPol1.
- indexA = indexElVertInPol2;
- indexB = indexElVertInPol1;
- } else if (fabs((*intPoints)[0][indexElVertInPol2]) <= 1.e-15) {
+ indexA = indexElVertInPol2;
+ indexB = indexElVertInPol1;
+ } else if (fabs((*intPoints)[0][indexElVertInPol2]) <= 1.e-15) {
- // (*intPoints)[0] lies in the face opposite to vertex
- // /ref indexElVertInPol2.
+ // (*intPoints)[0] lies in the face opposite to vertex
+ // /ref indexElVertInPol2.
- indexA = indexElVertInPol1;
- indexB = indexElVertInPol2;
- } else {
- ERROR_EXIT("couldn't determine A and B\n");
- }
+ indexA = indexElVertInPol1;
+ indexB = indexElVertInPol2;
+ } else {
+ ERROR_EXIT("couldn't determine A and B\n");
+ }
- /**
- * Sort the intersection points.
- */
+ /**
+ * Sort the intersection points.
+ */
- // (*intPoints)[0] is a neighbour of A (A has been constructed this way).
- indexS_0 = 0;
+ // (*intPoints)[0] is a neighbour of A (A has been constructed this way).
+ indexS_0 = 0;
- if (fabs((*intPoints)[1][indexB]) <= 1.e-15) {
+ if (fabs((*intPoints)[1][indexB]) <= 1.e-15) {
- // (*intPoints)[1] lies in the face opposite to B, thus is a neighbour
- // of A.
+ // (*intPoints)[1] lies in the face opposite to B, thus is a neighbour
+ // of A.
- indexS_1 = 1;
+ indexS_1 = 1;
- indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
+ indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
- if (fabs((*intPoints)[2][indexSecondFaceIntPoint0]) <= 1.e-15) {
+ if (fabs((*intPoints)[2][indexSecondFaceIntPoint0]) <= 1.e-15) {
- // (*intPoints)[2] is neighbour of (*intPoints)[0]
+ // (*intPoints)[2] is neighbour of (*intPoints)[0]
- indexS_2 = 3;
- indexS_3 = 2;
- } else {
+ indexS_2 = 3;
+ indexS_3 = 2;
+ } else {
- // (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
- // plane
+ // (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
+ // plane
- indexS_2 = 2;
- indexS_3 = 3;
- }
- } else if (fabs((*intPoints)[1][indexA]) <= 1.e-15) {
+ indexS_2 = 2;
+ indexS_3 = 3;
+ }
+ } else if (fabs((*intPoints)[1][indexA]) <= 1.e-15) {
- // (*intPoints)[1] lies in the face opposite to A
+ // (*intPoints)[1] lies in the face opposite to A
- indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
+ indexSecondFaceIntPoint0 = getIndexSecondFaceIntPoint0(indexB, dim);
- if (fabs((*intPoints)[1][indexSecondFaceIntPoint0]) <= 1.e-15) {
+ if (fabs((*intPoints)[1][indexSecondFaceIntPoint0]) <= 1.e-15) {
- // (*intPoints)[1] is neighbour of (*intPoints)[0], but isn't
- // neighbour of A
+ // (*intPoints)[1] is neighbour of (*intPoints)[0], but isn't
+ // neighbour of A
- indexS_3 = 1;
+ indexS_3 = 1;
- if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
+ if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
- // (*intPoints)[2] is neighbour of (*intPoints)[0] and neighbour of A
+ // (*intPoints)[2] is neighbour of (*intPoints)[0] and neighbour of A
- indexS_1 = 2;
- indexS_2 = 3;
- } else {
+ indexS_1 = 2;
+ indexS_2 = 3;
+ } else {
- // (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
- // plane
+ // (*intPoints)[2] is opposite to (*intPoints)[0] in the intersection
+ // plane
- indexS_1 = 3;
- indexS_2 = 2;
- }
- } else {
+ indexS_1 = 3;
+ indexS_2 = 2;
+ }
+ } else {
- // (*intPoints)[1] isn't neighbour of (*intPoints)[0], thus lies opposite
- // to (*intPoints)[0] in the intersection plane
+ // (*intPoints)[1] isn't neighbour of (*intPoints)[0], thus lies opposite
+ // to (*intPoints)[0] in the intersection plane
- indexS_2 = 1;
+ indexS_2 = 1;
- if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
+ if (fabs((*intPoints)[2][indexB]) <= 1.e-15) {
- // (*intPoints)[2] is neighbour of A
+ // (*intPoints)[2] is neighbour of A
- indexS_1 = 2;
- indexS_3 = 3;
- } else {
+ indexS_1 = 2;
+ indexS_3 = 3;
+ } else {
- // (*intPoints)[2] isn't neighbour of A
+ // (*intPoints)[2] isn't neighbour of A
- indexS_1 = 3;
- indexS_3 = 2;
+ indexS_1 = 3;
+ indexS_3 = 2;
+ }
}
+ } else {
+ ERROR_EXIT("IntPoint1 isn't either part of the face opposite to A nor part of the face opposite to B\n");
}
- } else {
- ERROR_EXIT("IntPoint1 isn't either part of the face opposite to A nor part of the face opposite to B\n");
+
+ /**
+ * For each subelement: Collect the vertices of the subelement in
+ * subElVertices, create a new SubElInfo for the subelement and
+ * store it in subElements.
+ *
+ * Note: The routines CompositeFEMOperator::getElementMatrix and
+ * CompositeFEMOperator::getElementVector expect the first vertice
+ * of a subelement to be a vertice of the corresponding element and
+ * not to be an intersection point.
+ */
+
+ // Create vertex A and vertex B.
+ vertexA[indexA] = 1.0;
+ vertexB[indexB] = 1.0;
+
+ // Subelement 1: A - B - S_0 - S_1
+ (*subElVertices)[0] = vertexA;
+ (*subElVertices)[1] = vertexB;
+ (*subElVertices)[2] = (*intPoints)[indexS_0];
+ (*subElVertices)[3] = (*intPoints)[indexS_1];
+ subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
+
+ // Subelement 2: B - S_0 - S_1 - S_2
+ (*subElVertices)[0] = vertexB;
+ (*subElVertices)[1] = (*intPoints)[indexS_0];
+ (*subElVertices)[2] = (*intPoints)[indexS_1];
+ (*subElVertices)[3] = (*intPoints)[indexS_2];
+ subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
+
+ // Subelement 3: B - S_0 - S_2 - S_3
+ (*subElVertices)[0] = vertexB;
+ (*subElVertices)[1] = (*intPoints)[indexS_0];
+ (*subElVertices)[2] = (*intPoints)[indexS_2];
+ (*subElVertices)[3] = (*intPoints)[indexS_3];
+ subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
+
+ TEST_EXIT( subElements.size() == 3 )("error in creating subelements");
+ numSubElements = 3;
+
+ DELETE subElVertices;
}
- /**
- * For each subelement: Collect the vertices of the subelement in
- * subElVertices, create a new SubElInfo for the subelement and
- * store it in subElements.
- *
- * Note: The routines CompositeFEMOperator::getElementMatrix and
- * CompositeFEMOperator::getElementVector expect the first vertice
- * of a subelement to be a vertice of the corresponding element and
- * not to be an intersection point.
- */
-
- // Create vertex A and vertex B.
- vertexA[indexA] = 1.0;
- vertexB[indexB] = 1.0;
-
- // Subelement 1: A - B - S_0 - S_1
- (*subElVertices)[0] = vertexA;
- (*subElVertices)[1] = vertexB;
- (*subElVertices)[2] = (*intPoints)[indexS_0];
- (*subElVertices)[3] = (*intPoints)[indexS_1];
- subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
-
- // Subelement 2: B - S_0 - S_1 - S_2
- (*subElVertices)[0] = vertexB;
- (*subElVertices)[1] = (*intPoints)[indexS_0];
- (*subElVertices)[2] = (*intPoints)[indexS_1];
- (*subElVertices)[3] = (*intPoints)[indexS_2];
- subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
-
- // Subelement 3: B - S_0 - S_2 - S_3
- (*subElVertices)[0] = vertexB;
- (*subElVertices)[1] = (*intPoints)[indexS_0];
- (*subElVertices)[2] = (*intPoints)[indexS_2];
- (*subElVertices)[3] = (*intPoints)[indexS_3];
- subElements.push_back( NEW SubElInfo(subElVertices, elInfo) );
-
- TEST_EXIT( subElements.size() == 3 )("error in creating subelements");
- numSubElements = 3;
-
- DELETE subElVertices;
}
diff --git a/AMDiS/compositeFEM/src/SubPolytope.h b/AMDiS/compositeFEM/src/SubPolytope.h
index d91a02028bf3869936b7aea60f482d76a11e8250..1b59a5888e2f4493fdfbda25630c24293ef637e0 100644
--- a/AMDiS/compositeFEM/src/SubPolytope.h
+++ b/AMDiS/compositeFEM/src/SubPolytope.h
@@ -8,195 +8,190 @@
#include "MemoryManager.h"
#include "Mesh.h"
#include "FixVec.h"
-
#include "SubElInfo.h"
-using namespace AMDiS;
-using namespace std;
-
-// ===========================================================================
-// === class SubPolytope =====================================================
-// ===========================================================================
-//
-// Class description:
-// The class SubPolytope holds the functionality for the division of a
-// subpolytope in subelements and contains a list of these subelements
-// in subElements. The number of subelements is given in numSubElements.
-//
-// The subpolytope of the element elInfo is initially given by the
-// intersection points (barycentric coordinates with respect to element)
-// resulting from the intersection of a "plane" with the element.
-// There is/are
-// 1 intersection point in 1-dimensional finite element space,
-// 2 intersection points in 2-dimensional finite element space,
-// 3 or 4 intersection points in 3-dimensional finite element space.
-// A detailed description of how the division into subelements works can be
-// found in the corresponding routines.
-//
-// Note: This functionality is restricted to the finite element dimensions
-// 1, 2 and 3.
-//
-// Note: The intersection points are expected to lie exactly on edges.
-// That means in the barycentric coordinate vector there are exactly
-// as many components equal to zero as are expected for points on edges.
-//
-// Main routines:
-// SubPolytope() - Creates the subelements for a subpolytope given by the
-// intersection points intPoints and stores them in
-// subElements.
-// createSubElementPolytopeIsSubElement1D()
-// - Creates the subelement in 1-dimensional finite element
-// space.
-// createSubElementPolytopeIsSubElement2D3D()
-// - Creates the subelement in 2-dimensional finite element
-// space and in 3-dimensional finite element space if there
-// are 3 intersection points.
-// createSubElementsForPolytope3D()
-// - Creates the 3 subelements in 3-dimensional finite element
-// space if there are 4 intersection points.
-// ============================================================================
-
-class SubPolytope
-{
- public:
- MEMORY_MANAGED(SubPolytope);
-
- /**
- * Constructor
- *
- * Divides the polytope produced by the intersection into subelements.
- * In dimensions 1 and 3 (case "4 intersection points") indexElVertInPol_
- * indicates which subpolytope to divide. The element vertice with index
- * indexElVertInPol_ is a vertice of the choosen subpolytope.
- * Here the standard vertice numeration is used (e.g. in dimension 3:
- * (1,0,0,0) - index 0, (0,1,0,0) - index 1, ...)
- * The subelements are stored in subElements.
- */
- SubPolytope(const ElInfo *elInfo_,
- VectorOfFixVecs<DimVec<double> > *intPoints_,
- int numIntPoints_,
- const int &indexElVertInPol_=0);
-
- /**
- * Destructor
- */
- ~SubPolytope()
+namespace AMDiS {
+
+ // ===========================================================================
+ // === class SubPolytope =====================================================
+ // ===========================================================================
+ //
+ // Class description:
+ // The class SubPolytope holds the functionality for the division of a
+ // subpolytope in subelements and contains a list of these subelements
+ // in subElements. The number of subelements is given in numSubElements.
+ //
+ // The subpolytope of the element elInfo is initially given by the
+ // intersection points (barycentric coordinates with respect to element)
+ // resulting from the intersection of a "plane" with the element.
+ // There is/are
+ // 1 intersection point in 1-dimensional finite element space,
+ // 2 intersection points in 2-dimensional finite element space,
+ // 3 or 4 intersection points in 3-dimensional finite element space.
+ // A detailed description of how the division into subelements works can be
+ // found in the corresponding routines.
+ //
+ // Note: This functionality is restricted to the finite element dimensions
+ // 1, 2 and 3.
+ //
+ // Note: The intersection points are expected to lie exactly on edges.
+ // That means in the barycentric coordinate vector there are exactly
+ // as many components equal to zero as are expected for points on edges.
+ //
+ // Main routines:
+ // SubPolytope() - Creates the subelements for a subpolytope given by the
+ // intersection points intPoints and stores them in
+ // subElements.
+ // createSubElementPolytopeIsSubElement1D()
+ // - Creates the subelement in 1-dimensional finite element
+ // space.
+ // createSubElementPolytopeIsSubElement2D3D()
+ // - Creates the subelement in 2-dimensional finite element
+ // space and in 3-dimensional finite element space if there
+ // are 3 intersection points.
+ // createSubElementsForPolytope3D()
+ // - Creates the 3 subelements in 3-dimensional finite element
+ // space if there are 4 intersection points.
+ // ============================================================================
+
+ class SubPolytope
{
- int i;
-
- for (i = 0; i < numSubElements; i++) {
- DELETE subElements[i];
+ public:
+ MEMORY_MANAGED(SubPolytope);
+
+ /**
+ * Constructor
+ *
+ * Divides the polytope produced by the intersection into subelements.
+ * In dimensions 1 and 3 (case "4 intersection points") indexElVertInPol_
+ * indicates which subpolytope to divide. The element vertice with index
+ * indexElVertInPol_ is a vertice of the choosen subpolytope.
+ * Here the standard vertice numeration is used (e.g. in dimension 3:
+ * (1,0,0,0) - index 0, (0,1,0,0) - index 1, ...)
+ * The subelements are stored in subElements.
+ */
+ SubPolytope(const ElInfo *elInfo_,
+ VectorOfFixVecs<DimVec<double> > *intPoints_,
+ int numIntPoints_,
+ const int &indexElVertInPol_ = 0);
+
+ /**
+ * Destructor
+ */
+ ~SubPolytope() {
+ for (int i = 0; i < numSubElements; i++) {
+ DELETE subElements[i];
+ }
+ };
+
+ /**
+ * Returns begin of vector subElements.
+ */
+ inline std::vector<SubElInfo *>::iterator getSubElementsBegin() {
+ return subElements.begin();
}
- };
-
- /**
- * Returns begin of vector subElements.
- */
- inline std::vector<SubElInfo *>::iterator getSubElementsBegin()
- {
- return subElements.begin();
- }
- /**
- * Returns end of vector subElements.
- */
- inline std::vector<SubElInfo *>::iterator getSubElementsEnd()
- {
- return subElements.end();
- }
+ /**
+ * Returns end of vector subElements.
+ */
+ inline std::vector<SubElInfo *>::iterator getSubElementsEnd() {
+ return subElements.end();
+ }
- /**
- * Returns num-th subelement of polytope.
- */
- inline SubElInfo* getSubElement(int num)
- {
- FUNCNAME("SubPolytope::getSubElement");
+ /**
+ * Returns num-th subelement of polytope.
+ */
+ inline SubElInfo* getSubElement(int num) {
+ FUNCNAME("SubPolytope::getSubElement()");
- TEST_EXIT(num <= numSubElements)("invalid index for subelement");
- return subElements[num];
- }
-
- /**
- * Returns the number of subelements of the polytope.
- */
- inline int getNumSubElements() { return numSubElements; }
-
- protected:
-
- /**
- * Checks whether the elements of intPoints are really intersection points,
- * i.e. lie on edges.
- */
- bool checkIntPoints();
-
- /**
- * In 1-dimensional space the intersection of an element always produces
- * two subelements. indexElVertInPol indicates which subelement to take.
- * The resulting subelement is created by this routine.
- */
- void createSubElementPolytopeIsSubElement1D(int indexElVertInPol);
-
- /**
- * In 2-dimensional space the intersection of an element always produces
- * a subelement. The same in 3-dimensional space, if the intersection has
- * three intersection points.
- * The resulting subelement is created by this routine.
- */
- void createSubElementPolytopeIsSubElement2D3D();
-
- /**
- * Routine used in createSubElementsForPolytope().
- *
- * The intersection point /ref intPoints[0] lies on an edge of element
- * and isn't a vertex of element. Thus it lies in exactly two faces.
- * In barycentric coordinates, lying in the face opposite to e.g. (0,1,0,0)
- * means that the second barycentric coordinate is equal to zero. We
- * give this face the index 1 (we start indexing by 0).
- * In this routine we know already one face containing intPoints[0],
- * namely indexFirstFace. The task of the routine is to get the
- * second face.
- */
- int getIndexSecondFaceIntPoint0(int indexFirstFace, int dim);
-
- /**
- * If in 3-dimensional space the intersection of an element has four
- * intersection points, the resulting polytope is no subelement, but it
- * can be divided into three subelements. This is done by this routine.
- */
- void createSubElementsForPolytope3D(int indexElVertInPol);
-
- protected:
-
- /**
- * elInfo of the element containing the polytope
- */
- const ElInfo *elInfo;
-
- /**
- * Intersection points with the element in barycentric coordinates with
- * respect to element
- */
- VectorOfFixVecs<DimVec<double> > *intPoints;
-
- /**
- * Number of intersection points
- */
- int numIntPoints;
-
- /**
- * List of the subelements of subpolytope
- */
- vector<SubElInfo *> subElements;
-
- /**
- * Number of subelements
- */
- int numSubElements;
-
- /**
- * Dimension of the polytope
- */
- int dim;
-};
+ TEST_EXIT(num <= numSubElements)("invalid index for subelement");
+ return subElements[num];
+ }
+
+ /**
+ * Returns the number of subelements of the polytope.
+ */
+ inline int getNumSubElements() {
+ return numSubElements;
+ }
+
+ protected:
+
+ /**
+ * Checks whether the elements of intPoints are really intersection points,
+ * i.e. lie on edges.
+ */
+ bool checkIntPoints();
+
+ /**
+ * In 1-dimensional space the intersection of an element always produces
+ * two subelements. indexElVertInPol indicates which subelement to take.
+ * The resulting subelement is created by this routine.
+ */
+ void createSubElementPolytopeIsSubElement1D(int indexElVertInPol);
+
+ /**
+ * In 2-dimensional space the intersection of an element always produces
+ * a subelement. The same in 3-dimensional space, if the intersection has
+ * three intersection points.
+ * The resulting subelement is created by this routine.
+ */
+ void createSubElementPolytopeIsSubElement2D3D();
+
+ /**
+ * Routine used in createSubElementsForPolytope().
+ *
+ * The intersection point /ref intPoints[0] lies on an edge of element
+ * and isn't a vertex of element. Thus it lies in exactly two faces.
+ * In barycentric coordinates, lying in the face opposite to e.g. (0,1,0,0)
+ * means that the second barycentric coordinate is equal to zero. We
+ * give this face the index 1 (we start indexing by 0).
+ * In this routine we know already one face containing intPoints[0],
+ * namely indexFirstFace. The task of the routine is to get the
+ * second face.
+ */
+ int getIndexSecondFaceIntPoint0(int indexFirstFace, int dim);
+
+ /**
+ * If in 3-dimensional space the intersection of an element has four
+ * intersection points, the resulting polytope is no subelement, but it
+ * can be divided into three subelements. This is done by this routine.
+ */
+ void createSubElementsForPolytope3D(int indexElVertInPol);
+
+ protected:
+
+ /**
+ * elInfo of the element containing the polytope
+ */
+ const ElInfo *elInfo;
+
+ /**
+ * Intersection points with the element in barycentric coordinates with
+ * respect to element
+ */
+ VectorOfFixVecs<DimVec<double> > *intPoints;
+
+ /**
+ * Number of intersection points
+ */
+ int numIntPoints;
+
+ /**
+ * List of the subelements of subpolytope
+ */
+ std::vector<SubElInfo *> subElements;
+
+ /**
+ * Number of subelements
+ */
+ int numSubElements;
+
+ /**
+ * Dimension of the polytope
+ */
+ int dim;
+ };
+}
#endif // AMDIS_SUBPOLYTOPE_H