From c8200a8f66f193e0ec0d983ff7a4ac1a61ab88bc Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Mon, 7 Dec 2009 13:03:33 +0000
Subject: [PATCH] - setting pointers as arguments to limit copying - use
fullmatrix to store levelset values at nodes - fix compile + clean up
---
Geo/MElementCut.cpp | 250 ++-
Solver/multiscaleLaplace.cpp | 2 +-
contrib/DiscreteIntegration/Integration3D.cpp | 1649 +++++++++--------
contrib/DiscreteIntegration/Integration3D.h | 316 ++--
contrib/DiscreteIntegration/recurCut.cpp | 530 ++----
contrib/DiscreteIntegration/recurCut.h | 12 +-
6 files changed, 1295 insertions(+), 1464 deletions(-)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index ac50ddb52e..667d49dae2 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -21,7 +21,7 @@ void MPolyhedron::_init()
if(_parts.size() == 0) return;
for(unsigned int i = 0; i < _parts.size(); i++) {
- if(_parts[i]->getVolume() * _parts[0]->getVolume() < 0)
+ if(_parts[i]->getVolume() * _parts[0]->getVolume() < 0.)
_parts[i]->revert();
for(int j = 0; j < 4; j++) {
int k;
@@ -217,47 +217,6 @@ void MPolygon::_initVertices()
}
}
- /*std::vector<MEdge>::iterator ite = edges.begin();
- MVertex *vINIT = ite->getVertex(0);
- _vertices.push_back(ite->getVertex(0));
- _vertices.push_back(ite->getVertex(1));
- edges.erase(ite);
-
- while(edges.size() > 1) {
- for(ite = edges.begin(); ite != edges.end(); ite++) {
- if(ite->getVertex(0) == _vertices.back()) {
- if(ite->getVertex(1) != vINIT) {
- _vertices.push_back(ite->getVertex(1));
- edges.erase(ite);
- }
- else {
- edges.erase(ite);
- ite = edges.begin();
- vINIT = ite->getVertex(0);
- _vertices.push_back(ite->getVertex(0));
- _vertices.push_back(ite->getVertex(1));
- edges.erase(ite);
- }
- break;
- }
- else if(ite->getVertex(1) == _vertices.back()) {
- if(ite->getVertex(0) != vINIT) {
- _vertices.push_back(ite->getVertex(0));
- edges.erase(ite);
- }
- else {
- edges.erase(ite);
- ite = edges.begin();
- vINIT = ite->getVertex(0);
- _vertices.push_back(ite->getVertex(0));
- _vertices.push_back(ite->getVertex(1));
- edges.erase(ite);
- }
- break;
- }
- }
- }*/
-
// innerVertices
for(unsigned int i = 0; i < _parts.size(); i++) {
for(int j = 0; j < 3; j++) {
@@ -427,14 +386,14 @@ void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
#if defined(HAVE_DINTEGRATION)
-static bool equalV(MVertex *v, const DI_Point &p)
+static bool equalV(MVertex *v, const DI_Point *p)
{
- return (fabs(v->x() - p.x()) < 1.e-15 &&
- fabs(v->y() - p.y()) < 1.e-15 &&
- fabs(v->z() - p.z()) < 1.e-15);
+ return (fabs(v->x() - p->x()) < 1.e-15 &&
+ fabs(v->y() - p->y()) < 1.e-15 &&
+ fabs(v->z() - p->z()) < 1.e-15);
}
-static int getElementVertexNum(DI_Point &p, MElement *e)
+static int getElementVertexNum(DI_Point *p, MElement *e)
{
for(int i = 0; i < e->getNumVertices(); i++)
if(equalV(e->getVertex(i), p))
@@ -442,7 +401,7 @@ static int getElementVertexNum(DI_Point &p, MElement *e)
return -1;
}
-static void assignPhysicals(GModel *GM, std::vector<int> gePhysicals, int reg, int dim,
+static void assignPhysicals(GModel *GM, std::vector<int> &gePhysicals, int reg, int dim,
std::map<int, std::map<int, std::string> > physicals[4])
{
for(unsigned int i = 0; i < gePhysicals.size(); i++){
@@ -495,8 +454,8 @@ static int getBorderTag(int lsTag, int count, int &maxTag, std::map<int, int> &b
typedef std::set<MVertex*,MVertexLessThanLexicographic> newVerticesContainer ;
-static void elementCutMesh(MElement *e, gLevelset *ls,
- std::map<int, std::map<int, double> > &verticesLs,
+static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
+ fullMatrix<double> &verticesLs,
GEntity *ge, GModel *GM, int &numEle, std::map<int, MVertex*> &vertexMap,
newVerticesContainer &newVertices,
std::map<int, std::vector<MElement*> > elements[10],
@@ -505,12 +464,12 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
std::map<int, int> newPhysTags[4],
std::map<int, int> borderElemTags[2],
std::map<int, int> borderPhysTags[2],
- std::vector<DI_CuttingPoint> &cp,
- std::vector<DI_Line> &lines,
- std::vector<DI_Triangle> &triangles,
- std::vector<DI_Quad> &quads,
- std::vector<DI_Tetra> &tetras,
- std::vector<DI_Hexa> &hexas)
+ DI_Point::Container &cp,
+ std::vector<DI_Line *> &lines,
+ std::vector<DI_Triangle *> &triangles,
+ std::vector<DI_Quad *> &quads,
+ std::vector<DI_Tetra *> &tetras,
+ std::vector<DI_Hexa *> &hexas)
{
int elementary = ge->tag();
int eType = e->getTypeForMSH();
@@ -518,8 +477,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
std::vector<int> gePhysicals = ge->physicals;
int integOrder = 1;
- std::vector<DI_IntegrationPoint> ipV;
- std::vector<DI_IntegrationPoint> ipS;
+ std::vector<DI_IntegrationPoint *> ipV;
+ std::vector<DI_IntegrationPoint *> ipS;
bool isCut = false;
unsigned int nbL = lines.size();
unsigned int nbTr = triangles.size();
@@ -560,6 +519,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
MElementFactory factory;
MElement *copy = factory.create(eType, vmv, ++numEle, ePart);
+ double **nodeLs = new double*[e->getNumPrimaryVertices()];
+ for(int i = 0; i < e->getNumPrimaryVertices(); i++) nodeLs[i] = new double[verticesLs.size2()];
+
switch (eType) {
case MSH_TET_4 :
case MSH_HEX_8 :
@@ -572,9 +534,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
- std::map<int, double> nodeLs[4];
- for(int i = 0; i < 4; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
- isCut = T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
}
else if(eType == MSH_HEX_8){
@@ -586,9 +547,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z(),
e->getVertex(6)->x(), e->getVertex(6)->y(), e->getVertex(6)->z(),
e->getVertex(7)->x(), e->getVertex(7)->y(), e->getVertex(7)->z());
- std::map<int, double> nodeLs[8];
- for(int i = 0; i < 8; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
- isCut = H.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
+ for(int i = 0; i < 8; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ isCut = H.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
hexas, tetras, quads, triangles, lines, 0, nodeLs);
}
else if(eType == MSH_PRI_6){
@@ -596,28 +556,26 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
- std::map<int, double> nodeLs[4];
- for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
- nodeLs[3] = verticesLs[e->getVertex(5)->getNum()];
- bool iC1 = T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ nodeLs[3] = &verticesLs(e->getVertex(5)->getIndex(),0);
+ bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
DI_Tetra T2(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z(),
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
- nodeLs[0] = verticesLs[e->getVertex(0)->getNum()];
- nodeLs[1] = verticesLs[e->getVertex(4)->getNum()];
- nodeLs[2] = verticesLs[e->getVertex(1)->getNum()];
- nodeLs[3] = verticesLs[e->getVertex(5)->getNum()];
- bool iC2 = T2.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ nodeLs[0] = &verticesLs(e->getVertex(0)->getIndex(),0);
+ nodeLs[1] = &verticesLs(e->getVertex(4)->getIndex(),0);
+ nodeLs[2] = &verticesLs(e->getVertex(1)->getIndex(),0);
+ nodeLs[3] = &verticesLs(e->getVertex(5)->getIndex(),0);
+ bool iC2 = T2.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
DI_Tetra T3(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z(),
e->getVertex(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
- nodeLs[0] = verticesLs[e->getVertex(0)->getNum()];
- for(int i = 1; i < 4; i++) nodeLs[i] = verticesLs[e->getVertex(i+2)->getNum()];
- bool iC3 = T3.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ for(int i = 1; i < 4; i++) nodeLs[i] = &verticesLs(e->getVertex(i+2)->getIndex(),0);
+ bool iC3 = T3.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
isCut = iC1 || iC2 || iC3;
}
@@ -626,31 +584,29 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
- std::map<int, double> nodeLs[4];
- for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
- nodeLs[3] = verticesLs[e->getVertex(4)->getNum()];
- bool iC1 = T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ nodeLs[3] = &verticesLs(e->getVertex(4)->getIndex(),0);
+ bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
DI_Tetra T2(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z(),
e->getVertex(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
- nodeLs[0] = verticesLs[e->getVertex(0)->getNum()];
- for(int i = 1; i < 4; i++) nodeLs[i] = verticesLs[e->getVertex(i+1)->getNum()];
- bool iC2 = T2.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ nodeLs[0] = &verticesLs(e->getVertex(0)->getIndex(),0);
+ for(int i = 1; i < 4; i++) nodeLs[i] = &verticesLs(e->getVertex(i+1)->getIndex(),0);
+ bool iC2 = T2.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
isCut = iC1 || iC2;
}
else if(eType == MSH_POLYH_){
- std::map<int, double> nodeLs[4];
for(int i = 0; i < e->getNumChildren(); i++) {
MTetrahedron *t = (MTetrahedron*) e->getChild(i);
DI_Tetra Tet(t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
t->getVertex(3)->x(), t->getVertex(3)->y(), t->getVertex(3)->z());
- for(int i = 0; i < 4; i++) nodeLs[i] = verticesLs[t->getVertex(i)->getNum()];
- bool iC = Tet.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(t->getVertex(i)->getIndex(),0);
+ bool iC = Tet.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
tetras, quads, triangles, 0, nodeLs);
isCut = (isCut || iC);
}
@@ -663,9 +619,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
for (unsigned int i = nbTe; i < tetras.size(); i++){
MVertex *mv[4] = {NULL, NULL, NULL, NULL};
for(int j = 0; j < 4; j++){
- int numV = getElementVertexNum(tetras[i].pt(j), e);
+ int numV = getElementVertexNum(tetras[i]->pt(j), e);
if (numV == -1){
- MVertex *newv = new MVertex(tetras[i].x(j), tetras[i].y(j), tetras[i].z(j), 0);
+ MVertex *newv = new MVertex(tetras[i]->x(j), tetras[i]->y(j), tetras[i]->z(j), 0);
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
if (!it.second) delete newv;
@@ -673,15 +629,15 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
if(it == vertexMap.end()) {
- mv[j] = new MVertex(tetras[i].x(j), tetras[i].y(j),
- tetras[i].z(j), 0, numV);
+ mv[j] = new MVertex(tetras[i]->x(j), tetras[i]->y(j),
+ tetras[i]->z(j), 0, numV);
vertexMap[numV] = mv[j];
}
else mv[j] = it->second;
}
}
MTetrahedron *mt = new MTetrahedron(mv[0], mv[1], mv[2], mv[3], ++numEle, ePart);
- if(tetras[i].lsTag() < 0)
+ if(tetras[i]->lsTag() < 0)
poly[0].push_back(mt);
else
poly[1].push_back(mt);
@@ -701,9 +657,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
}
else { // no cut
- int reg = getElementaryTag(tetras[nbTe].lsTag(), elementary, newElemTags[3]);
+ int reg = getElementaryTag(tetras[nbTe]->lsTag(), elementary, newElemTags[3]);
std::vector<int> phys;
- getPhysicalTag(tetras[nbTe].lsTag(), gePhysicals, phys, newPhysTags[3]);
+ getPhysicalTag(tetras[nbTe]->lsTag(), gePhysicals, phys, newPhysTags[3]);
if(eType == MSH_TET_4)
elements[4][reg].push_back(copy);
else if(eType == MSH_HEX_8)
@@ -720,9 +676,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
for (unsigned int i = nbTr; i < triangles.size(); i++){
MVertex *mv[3] = {NULL, NULL, NULL};
for(int j = 0; j < 3; j++){
- int numV = getElementVertexNum(triangles[i].pt(j), e);
+ int numV = getElementVertexNum(triangles[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(triangles[i].x(j), triangles[i].y(j), triangles[i].z(j), 0);
+ MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j), triangles[i]->z(j), 0);
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
if (!it.second) delete newv;
@@ -730,8 +686,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
if(it == vertexMap.end()) {
- mv[j] = new MVertex(triangles[i].x(j), triangles[i].y(j),
- triangles[i].z(j), 0, numV);
+ mv[j] = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
+ triangles[i]->z(j), 0, numV);
vertexMap[numV] = mv[j];
}
else mv[j] = it->second;
@@ -740,7 +696,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
MTriangle *tri;
if(p1 || p2) tri = new MTriangleBorder(mv[0], mv[1], mv[2], p1, p2, ++numEle, ePart);
else tri = new MTriangle(mv[0], mv[1], mv[2], ++numEle, ePart);
- int lsT = triangles[i].lsTag();
+ int lsT = triangles[i]->lsTag();
int c = elements[2].count(lsT) + elements[3].count(lsT);
// suppose that the surfaces have been cut before the volumes!
int reg = getBorderTag(lsT, c, newElemTags[2][0], borderElemTags[1]);
@@ -759,9 +715,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
DI_Triangle T(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z());
- std::map<int, double> nodeLs[3];
- for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
- isCut = T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
quads, triangles, lines, 0, nodeLs);
}
else if(eType == MSH_QUA_4){
@@ -769,20 +724,18 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z(),
e->getVertex(2)->x(), e->getVertex(2)->y(), e->getVertex(2)->z(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
- std::map<int, double> nodeLs[4];
- for(int i = 0; i < 4; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
- isCut = Q.cut(*ls, ipV, ipS, cp, integOrder,integOrder,integOrder,
+ for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ isCut = Q.cut(RPN, ipV, ipS, cp, integOrder,integOrder,integOrder,
quads, triangles, lines, 0, nodeLs);
}
else if(eType == MSH_POLYG_){
- std::map<int, double> nodeLs[3];
for(int i = 0; i < e->getNumChildren(); i++) {
MTriangle *t = (MTriangle*) e->getChild(i);
DI_Triangle Tri(t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z());
- for(int i = 0; i < 3; i++) nodeLs[i] = verticesLs[t->getVertex(i)->getNum()];
- bool iC = Tri.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(t->getVertex(i)->getIndex(),0);
+ bool iC = Tri.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
quads, triangles, lines, 0, nodeLs);
isCut = (isCut || iC);
}
@@ -795,9 +748,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
for (unsigned int i = nbTr; i < triangles.size(); i++){
MVertex *mv[3] = {NULL, NULL, NULL};
for(int j = 0; j < 3; j++){
- int numV = getElementVertexNum(triangles[i].pt(j), e);
+ int numV = getElementVertexNum(triangles[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(triangles[i].x(j), triangles[i].y(j), triangles[i].z(j), 0);
+ MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j), triangles[i]->z(j), 0);
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
if (!it.second) delete newv;
@@ -805,15 +758,15 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
if(it == vertexMap.end()) {
- mv[j] = new MVertex(triangles[i].x(j), triangles[i].y(j),
- triangles[i].z(j), 0, numV);
+ mv[j] = new MVertex(triangles[i]->x(j), triangles[i]->y(j),
+ triangles[i]->z(j), 0, numV);
vertexMap[numV] = mv[j];
}
else mv[j] = it->second;
}
}
MTriangle *mt = new MTriangle(mv[0], mv[1], mv[2], ++numEle, ePart);
- if(triangles[i].lsTag() < 0)
+ if(triangles[i]->lsTag() < 0)
poly[0].push_back(mt);
else
poly[1].push_back(mt);
@@ -833,9 +786,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
}
else { // no cut
- int reg = getElementaryTag(triangles[nbTr].lsTag(), elementary, newElemTags[2]);
+ int reg = getElementaryTag(triangles[nbTr]->lsTag(), elementary, newElemTags[2]);
std::vector<int> phys;
- getPhysicalTag(triangles[nbTr].lsTag(), gePhysicals, phys, newPhysTags[2]);
+ getPhysicalTag(triangles[nbTr]->lsTag(), gePhysicals, phys, newPhysTags[2]);
if(eType == MSH_TRI_3)
elements[2][reg].push_back(copy);
else if(eType == MSH_QUA_4)
@@ -848,9 +801,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
for (unsigned int i = nbL; i < lines.size(); i++){
MVertex *mv[2] = {NULL, NULL};
for(int j = 0; j < 2; j++){
- int numV = getElementVertexNum(lines[i].pt(j), e);
+ int numV = getElementVertexNum(lines[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(lines[i].x(j), lines[i].y(j), lines[i].z(j), 0);
+ MVertex *newv = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j), 0);
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
if (!it.second) delete newv;
@@ -858,8 +811,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
if(it == vertexMap.end()) {
- mv[j] = new MVertex(lines[i].x(j), lines[i].y(j),
- lines[i].z(j), 0, numV);
+ mv[j] = new MVertex(lines[i]->x(j), lines[i]->y(j),
+ lines[i]->z(j), 0, numV);
vertexMap[numV] = mv[j];
}
else mv[j] = it->second;
@@ -868,7 +821,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
MLine *lin;
if(p1 || p2) lin = new MLineBorder(mv[0], mv[1], p1, p2, ++numEle, ePart);
else lin = new MLine(mv[0], mv[1], ++numEle, ePart);
- int lsL = lines[i].lsTag();
+ int lsL = lines[i]->lsTag();
int c = elements[1].count(lsL);
// suppose that the lines have been cut before the surfaces!
int reg = getBorderTag(lsL, c, newElemTags[1][0], borderElemTags[0]);
@@ -883,17 +836,16 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
{
DI_Line L(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
e->getVertex(1)->x(), e->getVertex(1)->y(), e->getVertex(1)->z());
- std::map<int, double> nodeLs[2];
- for(int i = 0; i < 2; i++) nodeLs[i] = verticesLs[e->getVertex(i)->getNum()];
- isCut = L.cut(*ls, ipV, cp, integOrder, lines, 0, nodeLs);
+ for(int i = 0; i < 2; i++) nodeLs[i] = &verticesLs(e->getVertex(i)->getIndex(),0);
+ isCut = L.cut(RPN, ipV, cp, integOrder, lines, 0, nodeLs);
if(isCut) {
for (unsigned int i = nbL; i < lines.size(); i++){
MVertex *mv[2] = {NULL, NULL};
for(int j = 0; j < 2; j++){
- int numV = getElementVertexNum(lines[i].pt(j), e);
+ int numV = getElementVertexNum(lines[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(lines[i].x(j), lines[i].y(j), lines[i].z(j), 0);
+ MVertex *newv = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j), 0);
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
if (!it.second) delete newv;
@@ -901,24 +853,24 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
if(it == vertexMap.end()) {
- mv[j] = new MVertex(lines[i].x(j), lines[i].y(j), lines[i].z(j), 0, numV);
+ mv[j] = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j), 0, numV);
vertexMap[numV] = mv[j];
}
else mv[j] = it->second;
}
}
MLine *ml = new MLine(mv[0], mv[1], ++numEle, ePart);
- int reg = getElementaryTag(lines[i].lsTag(), elementary, newElemTags[1]);
+ int reg = getElementaryTag(lines[i]->lsTag(), elementary, newElemTags[1]);
std::vector<int> phys;
- getPhysicalTag(lines[i].lsTag(), gePhysicals, phys, newPhysTags[1]);
+ getPhysicalTag(lines[i]->lsTag(), gePhysicals, phys, newPhysTags[1]);
elements[1][reg].push_back(ml);
assignPhysicals(GM, phys, reg, 1, physicals);
}
}
else { // no cut
- int reg = getElementaryTag(lines[nbL].lsTag(), elementary, newElemTags[1]);
+ int reg = getElementaryTag(lines[nbL]->lsTag(), elementary, newElemTags[1]);
std::vector<int> phys;
- getPhysicalTag(lines[nbL].lsTag(), gePhysicals, phys, newPhysTags[1]);
+ getPhysicalTag(lines[nbL]->lsTag(), gePhysicals, phys, newPhysTags[1]);
elements[1][reg].push_back(copy);
assignPhysicals(GM, phys, reg, 1, physicals);
}
@@ -927,7 +879,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
case MSH_PNT :
{
DI_Point P(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z());
- P.computeLs(*ls);
+ P.computeLs(RPN.back());
int reg = getElementaryTag(P.lsTag(), elementary, newElemTags[0]);
std::vector<int> phys;
getPhysicalTag(P.lsTag(), gePhysicals, phys, newPhysTags[0]);
@@ -939,6 +891,14 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
Msg::Error("This type of element cannot be cut.");
return;
}
+
+ for(unsigned int i = 0; i < ipS.size(); i++)
+ delete ipS[i];
+ for(unsigned int i = 0; i < ipV.size(); i++)
+ delete ipV[i];
+ //for(int i = 0; i < e->getNumPrimaryVertices(); i++)
+ //if(nodeLs[i]) delete [] nodeLs[i];
+ delete [] nodeLs;
}
#endif
@@ -949,14 +909,14 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::map<int, std::map<int, std::string> > physicals[4])
{
#if defined(HAVE_DINTEGRATION)
+ int numVert = gm->indexMeshVertices(true);
+
GModel *cutGM = new GModel(gm->getName() + "_cut");
cutGM->setFileName(cutGM->getName());
newVerticesContainer newVertices;
- // std::set<MVertex*, MVertexLessThanLexicographic> newVertices;
std::vector<GEntity*> gmEntities;
-
gm->getEntities(gmEntities);
int numEle = gm->getNumMeshElements();
@@ -971,31 +931,39 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::vector<const gLevelset *> primitives;
ls->getPrimitives(primitives);
- std::map<int, std::map<int, double> > verticesLs;
+ fullMatrix<double> verticesLs(numVert + 1, primitives.size());
for(unsigned int i = 0; i < gmEntities.size(); i++) {
for(unsigned int j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
MVertex *vi = gmEntities[i]->getMeshVertex(j);
for(unsigned int k = 0; k < primitives.size(); k++) {
- verticesLs[vi->getNum()][primitives[k]->getTag()] = (*primitives[k])(vi->x(), vi->y(), vi->z());
+ verticesLs(vi->getIndex(), k) = (*primitives[k])(vi->x(), vi->y(), vi->z());
}
}
}
- std::vector<DI_CuttingPoint> cp;
- std::vector<DI_Line> lines;
- std::vector<DI_Triangle> triangles;
- std::vector<DI_Quad> quads;
- std::vector<DI_Tetra> tetras;
- std::vector<DI_Hexa> hexas;
+ DI_Point::Container cp;
+ std::vector<DI_Line *> lines;
+ std::vector<DI_Triangle *> triangles;
+ std::vector<DI_Quad *> quads;
+ std::vector<DI_Tetra *> tetras;
+ std::vector<DI_Hexa *> hexas;
+ std::vector<const gLevelset *> RPN;
+ ls->getRPN(RPN);
for(unsigned int i = 0; i < gmEntities.size(); i++) {
for(unsigned int j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
MElement *e = gmEntities[i]->getMeshElement(j);
e->setVolumePositive();
- elementCutMesh(e, ls, verticesLs, gmEntities[i], gm, numEle, vertexMap, newVertices,
+ elementCutMesh(e, RPN, verticesLs, gmEntities[i], gm, numEle, vertexMap, newVertices,
elements, physicals, newElemTags, newPhysTags, borderElemTags, borderPhysTags,
cp, lines, triangles, quads, tetras, hexas);
cutGM->getMeshPartitions().insert(e->getPartition());
}
+ for(DI_Point::Container::iterator it = cp.begin(); it != cp.end(); it++) delete *it;
+ for(unsigned int k = 0; k < lines.size(); k++) delete lines[k];
+ for(unsigned int k = 0; k < triangles.size(); k++) delete triangles[k];
+ for(unsigned int k = 0; k < quads.size(); k++) delete quads[k];
+ for(unsigned int k = 0; k < tetras.size(); k++) delete tetras[k];
+ for(unsigned int k = 0; k < hexas.size(); k++) delete hexas[k];
cp.clear(); lines.clear(); triangles.clear(); quads.clear(); tetras.clear(); hexas.clear();
}
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index 9324166fb3..b11d0812ae 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -11,7 +11,7 @@
#include "linearSystemCSR.h"
#ifdef HAVE_GMM
-#include "linearSystemGmm.h"
+#include "linearSystemGMM.h"
#endif
#include "linearSystemFull.h"
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 0ca3b33943..a38316cc18 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -26,8 +26,8 @@ inline double norm(const double *v) {
inline void vec (const double *x1, const double *x2, double *v) {
v[0] = x2[0] - x1[0]; v[1] = x2[1] - x1[1]; v[2] = x2[2] - x1[2];
}
-inline void vec (const DI_Point &p1, const DI_Point &p2, double *v) {
- v[0] = p2.x() - p1.x(); v[1] = p2.y() - p1.y(); v[2] = p2.z() - p1.z();
+inline void vec (const DI_Point *p1, const DI_Point *p2, double *v) {
+ v[0] = p2->x() - p1->x(); v[1] = p2->y() - p1->y(); v[2] = p2->z() - p1->z();
}
inline double sq2 (const double a) {return a * a;}
@@ -48,29 +48,29 @@ inline double det4 (double d11, double d12, double d13, double d14,
}
// distance
-inline double distance(const DI_Point &p1, const DI_Point &p2) {
- return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x())
- + (p1.y() - p2.y()) * (p1.y() - p2.y())
- + (p1.z() - p2.z()) * (p1.z() - p2.z()));
+inline double distance(const DI_Point *p1, const DI_Point *p2) {
+ return sqrt((p1->x() - p2->x()) * (p1->x() - p2->x())
+ + (p1->y() - p2->y()) * (p1->y() - p2->y())
+ + (p1->z() - p2->z()) * (p1->z() - p2->z()));
}
-inline double distance(const DI_CuttingPoint &p1, const DI_CuttingPoint &p2) {
- return sqrt((p1.x() - p2.x()) * (p1.x() - p2.x())
- + (p1.y() - p2.y()) * (p1.y() - p2.y())
- + (p1.z() - p2.z()) * (p1.z() - p2.z()));
+inline double distance(const DI_CuttingPoint *p1, const DI_CuttingPoint *p2) {
+ return sqrt((p1->x() - p2->x()) * (p1->x() - p2->x())
+ + (p1->y() - p2->y()) * (p1->y() - p2->y())
+ + (p1->z() - p2->z()) * (p1->z() - p2->z()));
}
// middle of 2 points
-inline DI_Point middle (const DI_Point &p1, const DI_Point &p2) {
- return DI_Point ((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2);
+inline DI_Point* middle (const DI_Point *p1, const DI_Point *p2) {
+ return new DI_Point ((p1->x() + p2->x()) / 2, (p1->y() + p2->y()) / 2, (p1->z() + p2->z()) / 2);
}
-inline DI_Point middle (const DI_Point &p1, const DI_Point &p2,
- const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
- return DI_Point ((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2, e, RPNi);
+inline DI_Point* middle (const DI_Point *p1, const DI_Point *p2,
+ const DI_Element *e, const std::vector<const gLevelset *> &RPNi) {
+ return new DI_Point ((p1->x() + p2->x()) / 2, (p1->y() + p2->y()) / 2, (p1->z() + p2->z()) / 2, e, RPNi);
}
// barycentre
-inline DI_Point barycenter (const DI_Element *el,
- const DI_Element *e, const std::vector<const gLevelset *> RPN) {
+inline DI_Point* barycenter (const DI_Element *el,
+ const DI_Element *e, const std::vector<const gLevelset *> &RPN) {
double x[3] = {0., 0., 0.};
int n;
for (n = 0; n < el->nbVert(); n++) {
@@ -78,7 +78,7 @@ inline DI_Point barycenter (const DI_Element *el,
x[1] += el->y(n);
x[2] += el->z(n);
}
- return DI_Point(x[0] / n,x[1] / n,x[2] / n, e, RPN);
+ return new DI_Point(x[0] / n, x[1] / n, x[2] / n, e, RPN);
}
// swap two integers
@@ -92,13 +92,13 @@ inline void swap(double &a, double &b) {
a = b; b = t;
}
// swap two points
-inline void swap (DI_Point &p1, DI_Point &p2) {
- DI_Point pt = p1;
- p1 = p2; p2 = pt;
+inline void swap (DI_Point **p1, DI_Point **p2) {
+ DI_Point *pt = *p1;
+ *p1 = *p2; *p2 = pt;
}
// return true if the 4 points are planar
-inline bool isPlanar (const DI_Point &p1, const DI_Point &p2, const DI_Point &p3, const DI_Point &p4) {
+inline bool isPlanar (const DI_Point *p1, const DI_Point *p2, const DI_Point *p3, const DI_Point *p4) {
double v12[3]; vec(p1, p2, v12);
double v13[3]; vec(p1, p3, v13);
double v14[3]; vec(p1, p4, v14);
@@ -108,14 +108,14 @@ inline bool isPlanar (const DI_Point &p1, const DI_Point &p2, const DI_Point &p3
return (c[0] == 0. && c[1] == 0. && c[2] == 0.);
}
// return true if the 3 points are linear
-inline bool isLinear (const DI_Point &p1, const DI_Point &p2, const DI_Point &p3) {
+inline bool isLinear (const DI_Point *p1, const DI_Point *p2, const DI_Point *p3) {
double v12[3]; vec(p1, p2, v12);
double v13[3]; vec(p1, p3, v13);
double c[3]; cross(v12, v13, c);
return (c[0] == 0. && c[1] == 0. && c[2] == 0.);
}
// return true if the 4 points form a quad in good ordering
-inline bool ordered4Nodes (const DI_Point &p1, const DI_Point &p2, const DI_Point &p3, const DI_Point &p4) {
+inline bool ordered4Nodes (const DI_Point *p1, const DI_Point *p2, const DI_Point *p3, const DI_Point *p4) {
double v21[3]; vec(p2, p1, v21);
double v23[3]; vec(p2, p3, v23);
double v24[3]; vec(p2, p4, v24);
@@ -130,52 +130,52 @@ inline bool ordered4Nodes (const DI_Point &p1, const DI_Point &p2, const DI_Poin
return true;
}
-bool isInQE (const DI_Triangle &t, const DI_QualError &QE) {
+bool isInQE (const DI_Triangle *t, const DI_QualError *QE) {
int b = 0;
for(int i = 0; i < 3; i++) {
for(int j = 0; j < 4; j++)
- if (t.pt(i).equal(QE.pt(j))) {b++; break;}
+ if (t->pt(i)->equal(QE->pt(j))) {b++; break;}
}
return (b == 3);
}
-bool belongsTo (const DI_Element &e, const DI_Element &E) {
+bool belongsTo (const DI_Element *e, const DI_Element *E) {
int b = 0;
- for(int k = 0; k < E.nbVert(); k++) {
- for(int l = 0; l < e.nbVert(); l++)
- if((e.pt(l)).equal(E.pt(k))) {b++; break;}
- if(b == e.nbVert()) return true;
+ for(int k = 0; k < E->nbVert(); k++) {
+ for(int l = 0; l < e->nbVert(); l++)
+ if((e->pt(l))->equal(E->pt(k))) {b++; break;}
+ if(b == e->nbVert()) return true;
}
return false;
}
-bool isLastLnInV (const std::vector<DI_Line> &v, const int i1 = 0) {
+bool isLastLnInV (const std::vector<DI_Line *> &v, const int i1 = 0) {
int b;
for (int i = i1; i < (int)v.size() - 1; i++) {
b = 0;
for (int k = 0; k < 2; k++)
for(int l = 0; l < 2; l++)
- if(v[i].pt(k).equal(v[v.size() - 1].pt(l))) {b++; break;}
+ if(v[i]->pt(k)->equal(v[v.size() - 1]->pt(l))) {b++; break;}
if(b == 2) return true;
}
return false;
}
-bool isLastTrInV (const std::vector<DI_Triangle> &v, const int i1 = 0) {
+bool isLastTrInV (const std::vector<DI_Triangle *> &v, const int i1 = 0) {
int b;
for (int i = i1; i < (int)v.size() - 1; i++) {
b = 0;
for (int k = 0; k < 3; k++)
for(int l = 0; l < 3; l++)
- if(v[i].pt(k).equal(v[v.size() - 1].pt(l))) {b++; break;}
+ if(v[i]->pt(k)->equal(v[v.size() - 1]->pt(l))) {b++; break;}
if(b == 3) return true;
}
return false;
}
-bool isLastQInV (const std::vector<DI_Quad> &v, const int i1 = 0) {
+bool isLastQInV (const std::vector<DI_Quad *> &v, const int i1 = 0) {
int b;
for (int i = i1; i < (int)v.size() - 1; i++) {
b = 0;
for(int k = 0; k < 4; k++)
for(int l = 0; l < 4; l++)
- if(v[i].pt(k).equal(v[v.size() - 1].pt(l))) {b++; break;}
+ if(v[i]->pt(k)->equal(v[v.size() - 1]->pt(l))) {b++; break;}
if(b == 4) return true;
}
return false;
@@ -192,10 +192,9 @@ int commonV (int &s11, int &s12, int &s21, int &s22) {
inline double adjustLs (double ls) {
return (fabs(ls) < ZERO_LS_TOL) ? 0. : ls;
}
-inline bool isCrossed (const DI_Point &p1, const DI_Point &p2) {
- double ls1 = p1.ls(), ls2 = p2.ls();
- if(ls1 * ls2 >= 0.) return false;
- return true;
+inline bool isCrossed (const DI_Point *p1, const DI_Point *p2) {
+ double ls1 = p1->ls(), ls2 = p2->ls();
+ return (ls1 * ls2 < 0.);
}
// return the index of the point with minimum x,y and z
@@ -220,7 +219,7 @@ int minimum(double *x, double *y, double *z, const int num) {
}
// Return the quality of a triangle
-double qualityTri(const DI_Point &p0, const DI_Point &p1, const DI_Point &p2) {
+double qualityTri(const DI_Point *p0, const DI_Point *p1, const DI_Point *p2) {
double a = distance(p0, p1);
double b = distance(p0, p2);
double c = distance(p1, p2);
@@ -231,14 +230,13 @@ double qualityTri(const DI_Point &p0, const DI_Point &p1, const DI_Point &p2) {
}
// Return the best quality for a quad cut into 2 triangles
-int bestQuality(const DI_Point &p1, const DI_Point &p2, const DI_Point &p3, const DI_Point &p4,
- DI_Triangle &t1, DI_Triangle &t2) {
+int bestQuality(const DI_Point *p1, const DI_Point *p2, const DI_Point *p3, const DI_Point *p4,
+ DI_Triangle **t1, DI_Triangle **t2) {
int cut = 0;
double qual11 = qualityTri(p1, p2, p3);
double qual12 = qualityTri(p1, p3, p4);
double qual21 = qualityTri(p1, p2, p4);
double qual22 = qualityTri(p2, p3, p4);
- //printf("quad, q11=%.18f,q12=%.18f,q21=%.18f,q22=%.18f\n",qual11,qual12,qual21,qual22);
double worst1 = std::min(qual11, qual12);
double worst2 = std::min(qual21, qual22);
if(worst1 - worst2 > EQUALITY_TOL)
@@ -253,22 +251,22 @@ int bestQuality(const DI_Point &p1, const DI_Point &p2, const DI_Point &p3, cons
else if(worst2 - worst1 > EQUALITY_TOL)
cut = 2;
else { // same min and max qualities
- double x[4] = {p1.x(), p2.x(), p3.x(), p4.x()};
- double y[4] = {p1.y(), p2.y(), p3.y(), p4.y()};
- double z[4] = {p1.z(), p2.z(), p3.z(), p4.z()};
+ double x[4] = {p1->x(), p2->x(), p3->x(), p4->x()};
+ double y[4] = {p1->y(), p2->y(), p3->y(), p4->y()};
+ double z[4] = {p1->z(), p2->z(), p3->z(), p4->z()};
int IND = minimum(x, y, z, 4);
if(IND == 0 || IND == 2) cut = 1;
else cut = 2;
}
}
if(cut == 1) {
- t1 = DI_Triangle(p1, p2, p3);
- t2 = DI_Triangle(p1, p3, p4);
+ *t1 = new DI_Triangle(p1, p2, p3);
+ *t2 = new DI_Triangle(p1, p3, p4);
return 1;
}
else {
- t1 = DI_Triangle(p1, p2, p4);
- t2 = DI_Triangle(p2, p3, p4);
+ *t1 = new DI_Triangle(p1, p2, p4);
+ *t2 = new DI_Triangle(p2, p3, p4);
return 2;
}
}
@@ -307,8 +305,8 @@ inline double qualityTet (double x1, double y1, double z1, double x2, double y2,
// Return the cutting of a pyramid cut into 2 tetras with the best quality faces
// (x1,x2,x3,x4 form the base in ccw and x5 is the summit)
-int bestQuality(const DI_Point &p1, const DI_Point &p2, const DI_Point &p3,
- const DI_Point &p4, const DI_Point &p5, DI_Tetra &t1, DI_Tetra &t2) {
+int bestQuality(const DI_Point *p1, const DI_Point *p2, const DI_Point *p3,
+ const DI_Point *p4, const DI_Point *p5, DI_Tetra **t1, DI_Tetra **t2) {
int cut = 0;
double qual11 = qualityTri(p1, p2, p3);
double qual12 = qualityTri(p1, p3, p4);
@@ -329,33 +327,33 @@ int bestQuality(const DI_Point &p1, const DI_Point &p2, const DI_Point &p3,
else if(worst2 - worst1 > EQUALITY_TOL)
cut = 2;
else { // same min and max qualities
- double x[4] = {p1.x(), p2.x(), p3.x(), p4.x()};
- double y[4] = {p1.y(), p2.y(), p3.y(), p4.y()};
- double z[4] = {p1.z(), p2.z(), p3.z(), p4.z()};
+ double x[4] = {p1->x(), p2->x(), p3->x(), p4->x()};
+ double y[4] = {p1->y(), p2->y(), p3->y(), p4->y()};
+ double z[4] = {p1->z(), p2->z(), p3->z(), p4->z()};
int IND = minimum(x, y, z, 4);
if(IND == 0 || IND == 2) cut = 1;
else cut = 2;
}
}
if(cut == 1) {
- t1 = DI_Tetra(p1, p2, p3, p5);
- t2 = DI_Tetra(p1, p3, p4, p5);
+ *t1 = new DI_Tetra(p1, p2, p3, p5);
+ *t2 = new DI_Tetra(p1, p3, p4, p5);
return 1;
}
else {
- t1 = DI_Tetra(p1, p2, p4, p5);
- t2 = DI_Tetra(p2, p3, p4, p5);
+ *t1 = new DI_Tetra(p1, p2, p4, p5);
+ *t2 = new DI_Tetra(p2, p3, p4, p5);
return 2;
}
}
// Return the cutting of a prism into 3 tetras with the best quality faces
-int bestQuality (const DI_Point &p0, const DI_Point &p1, const DI_Point &p2,
- const DI_Point &p3, const DI_Point &p4, const DI_Point &p5,
- DI_Tetra &t1, DI_Tetra &t2, DI_Tetra &t3, std::vector<DI_QualError> &QError) {
+int bestQuality (DI_Point *p0, DI_Point *p1, DI_Point *p2,
+ DI_Point *p3, DI_Point *p4, DI_Point *p5,
+ DI_Tetra **t1, DI_Tetra **t2, DI_Tetra **t3, std::vector<DI_QualError *> &QError) {
int cut3[3] = {0, 0, 0}; int cut = -1;
int fa[3][4] = {{0, 3, 4, 1}, {0, 2, 5, 3}, {1, 4, 5, 2}}; //faces
- DI_Point pt[6] = {p0, p1, p2, p3, p4, p5};
+ DI_Point* pt[6] = {p0, p1, p2, p3, p4, p5};
for(int f = 0; f < 3; f++){
int i1 = fa[f][0], i2 = fa[f][1], i3 = fa[f][2], i4 = fa[f][3];
double qual11 = qualityTri(pt[i1], pt[i2], pt[i3]);
@@ -377,9 +375,9 @@ int bestQuality (const DI_Point &p0, const DI_Point &p1, const DI_Point &p2,
else if(worst2 - worst1 > EQUALITY_TOL)
cut3[f] = 2;
else { // same min and max qualities
- double xc[4] = {pt[i1].x(), pt[i2].x(), pt[i3].x(), pt[i4].x()};
- double yc[4] = {pt[i1].y(), pt[i2].y(), pt[i3].y(), pt[i4].y()};
- double zc[4] = {pt[i1].z(), pt[i2].z(), pt[i3].z(), pt[i4].z()};
+ double xc[4] = {pt[i1]->x(), pt[i2]->x(), pt[i3]->x(), pt[i4]->x()};
+ double yc[4] = {pt[i1]->y(), pt[i2]->y(), pt[i3]->y(), pt[i4]->y()};
+ double zc[4] = {pt[i1]->z(), pt[i2]->z(), pt[i3]->z(), pt[i4]->z()};
int IND = minimum(xc, yc, zc, 4);
if(IND == 0 || IND == 2) cut3[f] = 1;
else cut3[f] = 2;
@@ -393,13 +391,13 @@ int bestQuality (const DI_Point &p0, const DI_Point &p1, const DI_Point &p2,
else cut = 2;
}
else{
- if(cut3[2] == 1) QError.push_back(DI_QualError(p1, p5, p2, p4));
+ if(cut3[2] == 1) QError.push_back(new DI_QualError(p1, p5, p2, p4));
cut = 3;
}
}
else {
if(cut3[1] == 1){
- if(cut3[2] == 2) QError.push_back(DI_QualError(p1, p5, p2, p4));
+ if(cut3[2] == 2) QError.push_back(new DI_QualError(p1, p5, p2, p4));
cut = 0;
}
else{
@@ -409,71 +407,72 @@ int bestQuality (const DI_Point &p0, const DI_Point &p1, const DI_Point &p2,
}
if(cut == 0) {
- t1 = DI_Tetra(p0, p1, p2, p5);
- t2 = DI_Tetra(p0, p1, p5, p3);
- t3 = DI_Tetra(p1, p5, p3, p4);
+ *t1 = new DI_Tetra(p0, p1, p2, p5);
+ *t2 = new DI_Tetra(p0, p1, p5, p3);
+ *t3 = new DI_Tetra(p1, p5, p3, p4);
return 1;
}
else if(cut == 1) {
- t1 = DI_Tetra(p0, p1, p2, p5);
- t2 = DI_Tetra(p0, p1, p5, p4);
- t3 = DI_Tetra(p0, p4, p5, p3);
+ *t1 = new DI_Tetra(p0, p1, p2, p5);
+ *t2 = new DI_Tetra(p0, p1, p5, p4);
+ *t3 = new DI_Tetra(p0, p4, p5, p3);
return 2;
}
else if(cut == 2) {
- t1 = DI_Tetra(p0, p1, p2, p4);
- t2 = DI_Tetra(p0, p4, p2, p5);
- t3 = DI_Tetra(p0, p4, p5, p3);
+ *t1 = new DI_Tetra(p0, p1, p2, p4);
+ *t2 = new DI_Tetra(p0, p4, p2, p5);
+ *t3 = new DI_Tetra(p0, p4, p5, p3);
return 3;
}
else if(cut == 3) {
- t1 = DI_Tetra(p0, p1, p2, p4);
- t2 = DI_Tetra(p0, p4, p2, p3);
- t3 = DI_Tetra(p2, p3, p4, p5);
+ *t1 = new DI_Tetra(p0, p1, p2, p4);
+ *t2 = new DI_Tetra(p0, p4, p2, p3);
+ *t3 = new DI_Tetra(p2, p3, p4, p5);
return 4;
}
else if(cut == 4) {
- t1 = DI_Tetra(p0, p1, p2, p3);
- t2 = DI_Tetra(p1, p2, p3, p4);
- t3 = DI_Tetra(p2, p3, p4, p5);
+ *t1 = new DI_Tetra(p0, p1, p2, p3);
+ *t2 = new DI_Tetra(p1, p2, p3, p4);
+ *t3 = new DI_Tetra(p2, p3, p4, p5);
return 5;
}
else if(cut == 5) {
- t1 = DI_Tetra(p0, p1, p2, p3);
- t2 = DI_Tetra(p1, p2, p3, p5);
- t3 = DI_Tetra(p1, p5, p3, p4);
+ *t1 = new DI_Tetra(p0, p1, p2, p3);
+ *t2 = new DI_Tetra(p1, p2, p3, p5);
+ *t3 = new DI_Tetra(p1, p5, p3, p4);
return 6;
}
}
// computes the intersection between a level set and a linear edge
-DI_Point Newton(const DI_Point &p1, const DI_Point &p2,
- const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
- double eps = -p1.ls() / (p2.ls() - p1.ls());
- DI_Point p(p1.x() + eps * (p2.x() - p1.x()), p1.y() + eps * (p2.y() - p1.y()), p1.z() + eps * (p2.z() - p1.z()));
- double pls = e->evalLs(p.x(), p.y(), p.z());
+DI_Point* Newton(const DI_Point *p1, const DI_Point *p2,
+ const DI_Element *e, const std::vector<const gLevelset *> &RPNi) {
+ double eps = -p1->ls() / (p2->ls() - p1->ls());
+ DI_Point* p = new DI_Point(p1->x() + eps * (p2->x() - p1->x()), p1->y() + eps * (p2->y() - p1->y()),
+ p1->z() + eps * (p2->z() - p1->z()));
+ double pls = e->evalLs(p->x(), p->y(), p->z());
// Newton iteration to find the intersection between the level set and the edge
while(fabs(pls) > EQUALITY_TOL){
- if(pls * p1.ls() < 0.) {
- eps = -pls / (p1.ls() - pls);
- p.move(p.x() + eps * (p1.x() - p.x()), p.y() + eps * (p1.y() - p.y()), p.z() + eps * (p1.z() - p.z()));
+ if(pls * p1->ls() < 0.) {
+ eps = -pls / (p1->ls() - pls);
+ p->move(p->x() + eps * (p1->x() - p->x()), p->y() + eps * (p1->y() - p->y()), p->z() + eps * (p1->z() - p->z()));
}
else {
- eps = -pls / (p2.ls() - pls);
- p.move(p.x() + eps * (p2.x() - p.x()), p.y() + eps * (p2.y() - p.y()), p.z() + eps * (p2.z() - p.z()));
+ eps = -pls / (p2->ls() - pls);
+ p->move(p->x() + eps * (p2->x() - p->x()), p->y() + eps * (p2->y() - p->y()), p->z() + eps * (p2->z() - p->z()));
}
- pls = e->evalLs(p.x(), p.y(), p.z());
+ pls = e->evalLs(p->x(), p->y(), p->z());
}
- p.computeLs(e, RPNi);
+ p->computeLs(e, RPNi);
return p;
}
// compute the position of the middle of a quadratic edge
//(intersection between the bisector of the linear edge and the levelset)
-DI_Point quadMidNode(const DI_Point &p1, const DI_Point &p2, const DI_Point &pf,
- const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
+DI_Point* quadMidNode(const DI_Point *p1, const DI_Point *p2, const DI_Point *pf,
+ const DI_Element *e, const std::vector<const gLevelset *> &RPNi) {
// middle of the edge between the 2 cutting points
- DI_Point midEN((p1.x() + p2.x()) / 2., (p1.y() + p2.y()) / 2., (p1.z() + p2.z()) / 2.);
+ DI_Point midEN((p1->x() + p2->x()) / 2., (p1->y() + p2->y()) / 2., (p1->z() + p2->z()) / 2.);
midEN.addLs(e);
// the bisector is computed in the plane of the face of the tetra [(x1,x2)x(x1,xf)]x(x1,x2)
double v12[3]; vec(p1, p2, v12);
@@ -492,15 +491,16 @@ DI_Point quadMidNode(const DI_Point &p1, const DI_Point &p2, const DI_Point &pf,
pt.move(midEN.x() - bisector[0], midEN.y() - bisector[1], midEN.z() - bisector[2]);
pt.changeLs(e->evalLs(pt.x(), pt.y(), pt.z()));
}
- return Newton(midEN, pt, e, RPNi);
+ DI_Point* qmn = Newton(&midEN, &pt, e, RPNi);
+ return qmn;
}
// ------------------------------------------------------------------------------------------------
// ------------------------------------------------------------------------------------------------
// DI_Point methods -------------------------------------------------------------------------------
-DI_Point::DI_Point (double x, double y, double z, gLevelset &ls) : x_(x), y_(y), z_(z) {
- addLs(ls(x, y, z));
+DI_Point::DI_Point (double x, double y, double z, gLevelset *ls) : x_(x), y_(y), z_(z) {
+ addLs((*ls)(x, y, z));
}
DI_Point & DI_Point::operator= (const DI_Point &rhs) {
if(this != &rhs){
@@ -521,7 +521,7 @@ void DI_Point::chooseLs (const gLevelset *Lsi) {
Ls.pop_back(); Ls.pop_back();
addLs(ls);
}
-void DI_Point::computeLs (const DI_Element *e, const std::vector<const gLevelset*> RPNi) {
+void DI_Point::computeLs (const DI_Element *e, const std::vector<const gLevelset*> &RPNi) {
int prim = 0; Ls.clear();
if(e->sizeLs() == 0) return;
for(int l = 0; l < (int)RPNi.size(); l++) {
@@ -531,16 +531,16 @@ void DI_Point::computeLs (const DI_Element *e, const std::vector<const gLevelset
chooseLs(RPNi[l]);
}
}
-void DI_Point::computeLs (const gLevelset &ls) {
+void DI_Point::computeLs (const gLevelset *ls) {
Ls.clear();
- addLs(ls(x_, y_, z_));
+ addLs((*ls)(x_, y_, z_));
}
-bool DI_Point::equal(const DI_Point &p) const {
- return (fabs(x() - p.x()) < EQUALITY_TOL && fabs(y() - p.y()) < EQUALITY_TOL && fabs(z() - p.z()) < EQUALITY_TOL);
+bool DI_Point::equal(const DI_Point *p) const {
+ return (fabs(x() - p->x()) < EQUALITY_TOL && fabs(y() - p->y()) < EQUALITY_TOL && fabs(z() - p->z()) < EQUALITY_TOL);
}
// DI_IntegrationPoint methods --------------------------------------------------------------------
-void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<const gLevelset*> RPN) {
+void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<const gLevelset*> &RPN) {
int prim = 0; std::vector<double> Ls;
for(int l = 0; l < (int)RPN.size(); l++) {
if(RPN[l]->isPrimitive())
@@ -556,26 +556,22 @@ void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<cons
}
// DI_CuttingPoint methods ------------------------------------------------------------------------
-DI_CuttingPoint::DI_CuttingPoint(const DI_Point &pt) : x_(pt.x()), y_(pt.y()), z_(pt.z()),
- xl_(pt.x()), yl_(pt.y()), zl_(pt.z()) {
- for(int i = 0; i < pt.sizeLs(); i++)
- addLs(pt.ls(i));
-}
-void DI_CuttingPoint::addLs (const double ls) {
- Ls.push_back(adjustLs(ls));
-}
-void DI_CuttingPoint::addLs (const DI_Element *e) {
- addLs(e->evalLs(x_, y_, z_));
+DI_CuttingPoint::DI_CuttingPoint(const DI_Point *pt) : DI_Point(pt->x(), pt->y(), pt->z()),
+ xl_(pt->x()), yl_(pt->y()), zl_(pt->z()) {
+ for(int i = 0; i < pt->sizeLs(); i++)
+ addLs(pt->ls(i));
}
-void DI_CuttingPoint::chooseLs (const gLevelset *Lsi) {
- if(Ls.size() < 2) return;
- double ls1 = Ls[Ls.size() - 2], ls2 = Ls[Ls.size() - 1];
- double ls = Lsi->choose(ls1, ls2);
- Ls.pop_back(); Ls.pop_back();
- addLs(ls);
-}
-bool DI_CuttingPoint::equal (const DI_CuttingPoint &p) const {
- return (fabs(x() - p.x()) < EQUALITY_TOL && fabs(y() - p.y()) < EQUALITY_TOL && fabs(z() - p.z()) < EQUALITY_TOL);
+
+// DI_PointLessThan methods -----------------------------------------------------------------------
+double DI_PointLessThan::tolerance = 1.e-6;
+bool DI_PointLessThan::operator()(const DI_Point *p1, const DI_Point *p2) const
+{
+ if(p1->x() - p2->x() > tolerance) return true;
+ if(p1->x() - p2->x() < -tolerance) return false;
+ if(p1->y() - p2->y() > tolerance) return true;
+ if(p1->y() - p2->y() < -tolerance) return false;
+ if(p1->z() - p2->z() > tolerance) return true;
+ return false;
}
// DI_Element methods -----------------------------------------------------------------------------
@@ -584,11 +580,11 @@ DI_Element::DI_Element(const DI_Element &cp) : lsTag_(cp.lsTag()), polOrder_(cp.
//printf("constructeur de copie d'element %d : ",cp.type()); cout << &cp << "->" << this << endl;
pts_ = new DI_Point* [cp.nbVert()];
for(int i = 0; i < cp.nbVert(); i++)
- pts_[i] = new DI_Point(cp.pt(i));
+ pts_[i] = new DI_Point(*cp.pt(i));
if(cp.getPolynomialOrder() > 1) {
mid_ = new DI_Point* [cp.nbMid()];
for(int i = 0; i < cp.nbMid(); i++)
- mid_[i] = new DI_Point(cp.mid(i));
+ mid_[i] = new DI_Point(*cp.mid(i));
}
else
mid_ = NULL;
@@ -610,7 +606,7 @@ DI_Element & DI_Element::operator= (const DI_Element &rhs){
if(this != &rhs) {
for(int i = 0; i < nbVert(); i++) {
delete pts_[i];
- pts_[i] = new DI_Point(rhs.pt(i));
+ pts_[i] = new DI_Point(*rhs.pt(i));
}
for(int i = 0; i < nbMid(); i++) {
delete mid_[i];
@@ -621,7 +617,7 @@ DI_Element & DI_Element::operator= (const DI_Element &rhs){
mid_ = new DI_Point*[rhs.nbMid()];
}
for(int i = 0; i < rhs.nbMid(); i++)
- mid_[i] = new DI_Point(rhs.mid(i));
+ mid_[i] = new DI_Point(*rhs.mid(i));
polOrder_ = rhs.getPolynomialOrder();
integral_ = rhs.integral();
lsTag_ = rhs.lsTag();
@@ -651,7 +647,7 @@ void DI_Element::setPolynomialOrder (int o) {
printf("Order %d line function space not implemented ", o); print();
}
}
-void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
+void DI_Element::setPolynomialOrder (int o, const DI_Element *e, const std::vector<const gLevelset *> &RPNi) {
if(polOrder_ == o) return;
for(int i = 0; i < nbMid(); i++)
delete mid_[i];
@@ -687,12 +683,12 @@ void DI_Element::addLs (const DI_Element *e) {
for(int i = 0; i < nbMid(); i++)
mid_[i]->addLs(e);
}
-void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls) {
+void DI_Element::addLs (const DI_Element *e, const gLevelset *Ls) {
if(type() != e->type()) {
printf("Error : addLs with element of different type\n");
}
for(int j = 0; j < nbVert(); ++j) {
- double ls = Ls(e->x(j), e->y(j), e->z(j));
+ double ls = (*Ls)(e->x(j), e->y(j), e->z(j));
pts_[j]->addLs(ls);
}
for(int j = 0; j < nbMid(); ++j) {
@@ -701,25 +697,24 @@ void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls) {
double xc = 0, yc = 0, zc = 0;
for(int k = 0; k < n; k++){
xc += e->x(s[k]); yc += e->y(s[k]); zc += e->z(s[k]); }
- double ls = Ls(xc/n, yc/n, zc/n);
+ double ls = (*Ls)(xc/n, yc/n, zc/n);
mid_[j]->addLs(ls);
}
}
-void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls, std::map<int, double> nodeLs[]) {
- int primTag = Ls.getTag();
- if(!nodeLs || !nodeLs[0].count(primTag)) {
+void DI_Element::addLs (const DI_Element *e, const gLevelset *Ls, int iLs, double **nodeLs) {
+ if(!nodeLs || !nodeLs[0][iLs]) {
addLs(e, Ls);
return;
}
for(int i = 0; i < nbVert(); i++)
- pts_[i]->addLs(nodeLs[i][primTag]);
+ pts_[i]->addLs(nodeLs[i][iLs]);
for(int j = 0; j < nbMid(); ++j) {
std::vector<int> s(nbVert()); int n;
e->midV(j, &s[0], n);
double xc = 0, yc = 0, zc = 0;
for(int k = 0; k < n; k++){
xc += e->x(s[k]); yc += e->y(s[k]); zc += e->z(s[k]); }
- double ls = Ls(xc/n, yc/n, zc/n);
+ double ls = (*Ls)(xc/n, yc/n, zc/n);
mid_[j]->addLs(ls);
}
}
@@ -738,12 +733,12 @@ void DI_Element::clearLs() {
mid_[i]->clearLs();
}
bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
- const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
+ const DI_Element *e, const std::vector<const gLevelset *> &RPNi) {
if(edge >= nbEdg()) {printf("wrong number (%d) for quadratic edge for a ", edge); print(); return false;}
int s1, s2; vert(edge, s1, s2);
int order0 = getPolynomialOrder();
if(order0 == 1) setPolynomialOrder(2, e, RPNi);
- double dist = distance(mid(edge), *xm);
+ double dist = distance(mid(edge), xm);
double sideLength = distance(pt(s1), pt(s2));
if (dist / sideLength < 1.e-5) {
if(order0 == 1) setPolynomialOrder(1);
@@ -763,9 +758,9 @@ bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
computeIntegral();
return true;
}
-bool DI_Element::contain (const DI_Point &p) const {
+bool DI_Element::contain (const DI_Point *p) const {
for(int i = 0; i < nbVert(); i++){
- if(p.equal(pt(i)))
+ if(p->equal(pt(i)))
return true;
}
switch(getDim()){
@@ -805,22 +800,19 @@ bool DI_Element::contain (const DI_Element *e) const {
if(!contain(e->pt(p))) return false;
return true;
}
-DI_Point DI_Element::mappingP (DI_Point &pt) const {
- double e[3]; evalC(pt.x(), pt.y(), pt.z(), e);
- pt.move(e[0], e[1], e[2]);
- return pt;
+void DI_Element::mappingP (DI_Point *pt) const {
+ double e[3]; evalC(pt->x(), pt->y(), pt->z(), e);
+ pt->move(e[0], e[1], e[2]);
}
-DI_IntegrationPoint DI_Element::mappingIP (DI_IntegrationPoint &in) const {
- double e[3]; evalC(in.x(), in.y(), in.z(), e);
- double w = in.weight() * integral() / refIntegral();
- in.move(e[0], e[1], e[2]);
- in.setWeight(w);
- return in;
+void DI_Element::mappingIP (DI_IntegrationPoint *in) const {
+ double e[3]; evalC(in->x(), in->y(), in->z(), e);
+ double w = in->weight() * integral() / refIntegral();
+ in->move(e[0], e[1], e[2]);
+ in->setWeight(w);
}
-DI_CuttingPoint DI_Element::mappingCP (DI_CuttingPoint &cp) const {
- double e[3]; evalC(cp.x(), cp.y(), cp.z(), e);
- cp.move(e[0], e[1], e[2]);
- return cp;
+void DI_Element::mappingCP (DI_CuttingPoint *cp) const {
+ double e[3]; evalC(cp->x(), cp->y(), cp->z(), e);
+ cp->move(e[0], e[1], e[2]);
}
void DI_Element::mappingEl (DI_Element *el) const {
double s[3];
@@ -835,36 +827,40 @@ void DI_Element::mappingEl (DI_Element *el) const {
el->computeIntegral();
}
void DI_Element::integrationPoints (const int polynomialOrder, const DI_Element *loc,
- const DI_Element *e, const std::vector<const gLevelset *> RPN,
- std::vector<DI_IntegrationPoint> &ip) const {
- std::vector<DI_IntegrationPoint> ip_ref;
+ const DI_Element *e, const std::vector<const gLevelset *> &RPN,
+ std::vector<DI_IntegrationPoint *> &ip) const {
+ std::vector<DI_IntegrationPoint *> ip_ref;
getRefIntegrationPoints(polynomialOrder, ip_ref);
for (int j = 0; j < (int)ip_ref.size(); j++) {
- DI_IntegrationPoint IPloc = ip_ref[j];
- loc->mappingIP(IPloc);
+ DI_IntegrationPoint IPloc (*ip_ref[j]);
+ loc->mappingIP(&IPloc);
mappingIP(ip_ref[j]);
- ip_ref[j].addLocC(IPloc.x(), IPloc.y(), IPloc.z());
- DI_IntegrationPoint pp = IPloc; pp.computeLs(e, RPN);
- //ip_ref[j].setLs(loc->evalLs(IPloc.x(), IPloc.y(), IPloc.z()));//levelset computed in the subelement FIXME check sign!
- //ip_ref[j].setLs(evalLs(ip_ref[j].x(), ip_ref[j].y(), ip_ref[j].z()));
- ip_ref[j].setLs(pp.ls());
- //printf("pt loc : "); IPloc.print(); printf("pt reel : "); ip_ref[j].print();
- //printf("el loc : "); loc->printls(); printf("el reel : "); printls(); printf("el e : "); e->printls();
- //printf("ls_loc=%g ls_reel=%g ls_e=%g\n\n", loc->evalLs(IPloc.x(), IPloc.y(), IPloc.z()), evalLs(ip_ref[j].x(), ip_ref[j].y(), ip_ref[j].z()), pp.ls());
+ ip_ref[j]->addLocC(IPloc.x(), IPloc.y(), IPloc.z());
+ DI_IntegrationPoint pp(IPloc); pp.computeLs(e, RPN);
+ //ip_ref[j]->setLs(loc->evalLs(IPloc.x(), IPloc.y(), IPloc.z()));//levelset computed in the subelement FIXME check sign!
+ //ip_ref[j]->setLs(evalLs(ip_ref[j]->x(), ip_ref[j]->y(), ip_ref[j]->z()));
+ ip_ref[j]->setLs(pp.ls());
+ /*printf("pt loc : "); IPloc.print(); printf("pt reel : "); ip_ref[j]->print();
+ printf("el loc : "); loc->printls(); printf("el reel : "); printls(); printf("el e : "); e->printls();
+ printf("ls_loc=%g ls_reel=%g ls_e=%g\n\n", loc->evalLs(IPloc.x(), IPloc.y(), IPloc.z()),
+ evalLs(ip_ref[j]->x(), ip_ref[j]->y(), ip_ref[j]->z()), pp.ls());*/
ip.push_back(ip_ref[j]);
}
}
-void DI_Element::getCuttingPoints (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
- std::vector<DI_CuttingPoint> &cp) const {
+void DI_Element::getCuttingPoints (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_CuttingPoint *> &cp) const {
int s1, s2;
for(int i = 0; i < nbEdg(); i++){
vert(i, s1, s2);
if (isCrossed(pt(s1), pt(s2))) {
- DI_Point p = Newton(pt(s1), pt(s2), e, RPNi);
- cp.push_back(DI_CuttingPoint(p));
+ DI_Point* p = Newton(pt(s1), pt(s2), e, RPNi);
+ cp.push_back(new DI_CuttingPoint(p));
+ delete p;
}
- if(ls(s1) == 0.)
- cp.push_back(DI_CuttingPoint(pt(s1)));
+ }
+ for(int i = 0; i < nbVert(); i++){
+ if(ls(i) == 0.)
+ cp.push_back(new DI_CuttingPoint(pt(i)));
}
}
void DI_Element::evalC (const double u, const double v, const double w, double *ev, int order) const {
@@ -938,34 +934,34 @@ double DI_Element::detJ (const double u, const double v, const double w) const {
}
}
// set the lsTag to +1 if the element is inside the domain (lsTag is -1 by default)
-void DI_Element::computeLsTagDom(const DI_Element *e, const std::vector<const gLevelset *> RPN) {
+void DI_Element::computeLsTagDom(const DI_Element *e, const std::vector<const gLevelset *> &RPN) {
for(int i = 0; i < nbVert(); i++) {
- if(pt(i).isOutsideDomain())
+ if(pt(i)->isOutsideDomain())
return;
- if(pt(i).isInsideDomain())
+ if(pt(i)->isInsideDomain())
{setLsTag(1); return;}
}
- DI_Point bar = barycenter(this, e, RPN);
- if(bar.isOutsideDomain())
- return;
- if(bar.isInsideDomain())
- {setLsTag(1); return;}
+ DI_Point* bar = barycenter(this, e, RPN);
+ if(bar->isOutsideDomain())
+ {delete bar; return;}
+ if(bar->isInsideDomain())
+ {setLsTag(1); delete bar; return;}
for(int i = 0; i < nbVert(); i++) {
- DI_Point mid = middle (pt(i), bar, e, RPN);
- if(mid.isOutsideDomain())
- return;
- if(mid.isInsideDomain())
- {setLsTag(1); return;}
+ DI_Point* mid = middle (pt(i), bar, e, RPN);
+ delete bar;
+ if(mid->isOutsideDomain())
+ {delete mid; return;}
+ if(mid->isInsideDomain())
+ {setLsTag(1); delete mid; return;}
}
printf("Error : Unable to determine the sign of the element : \n");
printf("Parent element : "); e->printls();
printf("Element : "); printls();
- return;
}
// set the lsTag to -1 if the element is not on the border of the domain
-void DI_Element::computeLsTagBound(std::vector<DI_Hexa> &hexas, std::vector<DI_Tetra> &tetras){
+void DI_Element::computeLsTagBound(std::vector<DI_Hexa *> &hexas, std::vector<DI_Tetra *> &tetras){
for(int i = 0; i < nbVert(); i++) {
- if(!pt(i).isOnBorder()) {
+ if(!pt(i)->isOnBorder()) {
setLsTag(-1);
return;
}
@@ -973,9 +969,9 @@ void DI_Element::computeLsTagBound(std::vector<DI_Hexa> &hexas, std::vector<DI_T
/*DI_Element *e1 = NULL, *e2 = NULL;
for(int i = 0; i < (int)tetras.size(); i++){
- if(belongsTo(*this, tetras[i])) {
- if(!e1) e1 = &tetras[i];
- else {e2 = &tetras[i]; break;}
+ if(belongsTo(this, tetras[i])) {
+ if(!e1) e1 = tetras[i];
+ else {e2 = tetras[i]; break;}
}
}
if(e1 && e2) {
@@ -983,14 +979,14 @@ void DI_Element::computeLsTagBound(std::vector<DI_Hexa> &hexas, std::vector<DI_T
return;
}
for(int i = 0; i < (int)hexas.size(); i++){
- if(belongsTo(*this, hexas[i])) {
- if(!e1) e1 = &hexas[i];
- else {e2 = &hexas[i]; break;}
+ if(belongsTo(this, hexas[i])) {
+ if(!e1) e1 = hexas[i];
+ else {e2 = hexas[i]; break;}
}
}
if(e1 && e2 && e1->lsTag() * e2->lsTag() > 0.) setLsTag(-1);*/
}
-void DI_Element::computeLsTagBound(std::vector<DI_Quad> &quads, std::vector<DI_Triangle> &triangles){
+void DI_Element::computeLsTagBound(std::vector<DI_Quad *> &quads, std::vector<DI_Triangle *> &triangles){
for(int i = 0; i < nbVert(); i++) {
if(ls(i) != 0.) {
setLsTag(-1);
@@ -1001,10 +997,10 @@ void DI_Element::computeLsTagBound(std::vector<DI_Quad> &quads, std::vector<DI_T
/*DI_Element *e1 = NULL, *e2 = NULL, *et = NULL;
int NT = triangles.size(), NQ = quads.size();
for(int i = 0; i < NT + NQ; i++){
- if((i < NT && belongsTo(*this, triangles[i])) ||
- (i >= NT && belongsTo(*this, quads[i-NT]))) {
- if(i < NT) et = &triangles[i];
- else et = &quads[i-NT];
+ if((i < NT && belongsTo(this, triangles[i])) ||
+ (i >= NT && belongsTo(this, quads[i-NT]))) {
+ if(i < NT) et = triangles[i];
+ else et = quads[i-NT];
if(!e1) e1 = et;
else {e2 = et; break;}
}
@@ -1044,6 +1040,19 @@ void DI_Element::printls() const {
printf("tag=%d\n", lsTag());
}
+double DI_ElementLessThan::tolerance = 1.e-6;
+bool DI_ElementLessThan::operator()(const DI_Element *e1, const DI_Element *e2) const
+{
+ for(int i = 0; i < e1->nbVert(); i++) {
+ if(e1->x(i) - e2->x(i) > tolerance) return true;
+ if(e1->x(i) - e2->x(i) < -tolerance) return false;
+ if(e1->y(i) - e2->y(i) > tolerance) return true;
+ if(e1->y(i) - e2->y(i) < -tolerance) return false;
+ if(e1->z(i) - e2->z(i) > tolerance) return true;
+ return false;
+ }
+}
+
// DI_Line methods --------------------------------------------------------------------------------
void DI_Line::computeIntegral() {
switch (getPolynomialOrder()) {
@@ -1416,60 +1425,62 @@ void DI_Hexa::getGradShapeFunctions (const double u, const double v, const doubl
// Split a reference line cut by a level set into sublines
void DI_Line::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cuttingPoints) const
+ std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
{
if (!isCrossed(pt(0), pt(1))) {
- subLines.push_back(*this);
+ subLines.push_back(new DI_Line(*this));
return;
}
// compute the intersection between the line and the level set
- DI_Point U = Newton(pt(0), pt(1), e, RPNi);
- cuttingPoints.push_back(DI_CuttingPoint(U));
+ DI_Point *U = Newton(pt(0), pt(1), e, RPNi);
+ cuttingPoints.push_back(new DI_CuttingPoint(U));
// line cut => split into 2 sublines
- subLines.push_back(DI_Line(pt(0), U, lsTag()));
- subLines.push_back(DI_Line(U, pt(1), lsTag()));
+ subLines.push_back(new DI_Line(pt(0), U, lsTag()));
+ subLines.push_back(new DI_Line(U, pt(1), lsTag()));
+ delete U;
}
// Split a triangle into 4 sub-triangles
-void DI_Triangle::splitIntoSubTriangles (std::vector<DI_Triangle> &subTriangles) const {
- DI_Point p01 = middle(pt(0), pt(1));
- DI_Point p02 = middle(pt(0), pt(2));
- DI_Point p12 = middle(pt(1), pt(2));
+void DI_Triangle::splitIntoSubTriangles (std::vector<DI_Triangle *> &subTriangles) const {
+ DI_Point *p01 = middle(pt(0), pt(1));
+ DI_Point *p02 = middle(pt(0), pt(2));
+ DI_Point *p12 = middle(pt(1), pt(2));
- subTriangles.push_back(DI_Triangle(p01, p02, p12)); // 01-02-12
- subTriangles.push_back(DI_Triangle(pt(0), p01, p02)); // 0-01-02
- subTriangles.push_back(DI_Triangle(pt(1), p01, p12)); // 1-01-12
- subTriangles.push_back(DI_Triangle(pt(2), p02, p12)); // 2-02-12
+ subTriangles.push_back(new DI_Triangle(p01, p02, p12)); // 01-02-12
+ subTriangles.push_back(new DI_Triangle(pt(0), p01, p02)); // 0-01-02
+ subTriangles.push_back(new DI_Triangle(pt(1), p01, p12)); // 1-01-12
+ subTriangles.push_back(new DI_Triangle(pt(2), p02, p12)); // 2-02-12
+ delete p01; delete p02; delete p12;
}
// Split a reference triangle cut by a level set into subtriangles
void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cuttingPoints) const
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
{
bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // if true, use quadratic sub triangles
int LStag = RPNi.back()->getTag();
int nbZe = 0, ze[2];
for(int i = 0; i < 3; i++)
- if(pt(i).isOnBorder())
+ if(pt(i)->isOnBorder())
ze[nbZe++] = i;
for(int i = 0; i < nbZe; i++)
- cuttingPoints.push_back(DI_CuttingPoint(pt(ze[i])));
+ cuttingPoints.push_back(new DI_CuttingPoint(pt(ze[i])));
// !! we add these points several times => remove later
if (!isCrossed(pt(0), pt(1)) && !isCrossed(pt(0), pt(2)) && !isCrossed(pt(1), pt(2))) {
- subTriangles.push_back(*this);
+ subTriangles.push_back(new DI_Triangle(*this));
if(nbZe == 2)
- surfLines.push_back(DI_Line(pt(ze[0]), pt(ze[1]), LStag));
+ surfLines.push_back(new DI_Line(pt(ze[0]), pt(ze[1]), LStag));
// !! we add these lines twice => remove the second later
return;
}
- DI_Point U[2]; // cutting points (max2)
- int IND[2]; // ids of edges cut
- int COUNT = 0; // number of edges cut
+ DI_Point *U[2]; // cutting points (max2)
+ int IND[2]; // ids of edges cut
+ int COUNT = 0; // number of edges cut
// compute the intersections between the edges of the triangle and the level set
for(int i = 0; i < 3; i++){
@@ -1481,7 +1492,7 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevel
}
for(int i = 0; i < COUNT; i++)
- cuttingPoints.push_back(DI_CuttingPoint(U[i]));
+ cuttingPoints.push_back(new DI_CuttingPoint(U[i]));
// Adjust the indices of the nodes in order to have the same pattern for each case with the same number of edges cut;
// compute the position of the middle nodes on the quadratic edges of the sub triangles;
@@ -1492,21 +1503,22 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevel
int s2 = (s1 + 1) % 3;
int s3 = (s2 + 1) % 3;
- DI_Triangle t1 = DI_Triangle(pt(s3), pt(s1), U[0], lsTag()); //t1.print();
- DI_Triangle t2 = DI_Triangle(pt(s2), pt(s3), U[0], lsTag()); //t2.print();
- DI_Line ln = DI_Line(U[0], pt(s3), LStag); //ln.print();
+ DI_Triangle *t1 = new DI_Triangle(pt(s3), pt(s1), U[0], lsTag());
+ DI_Triangle *t2 = new DI_Triangle(pt(s2), pt(s3), U[0], lsTag());
+ DI_Line *ln = new DI_Line(U[0], pt(s3), LStag);
if(quadT){
- DI_Point midEN2 = quadMidNode(U[0], pt(s3), pt(s1), e, RPNi); // intersection between ls and the bissector between the cutting points
- cuttingPoints.push_back(DI_CuttingPoint(midEN2));
- //midEN2.printls();
- t1.addQuadEdge (2, &midEN2, e, RPNi); //t1.printls();
- t2.addQuadEdge (1, &midEN2, e, RPNi); //t2.printls();
- ln.addQuadEdge (0, &midEN2, e, RPNi); //ln.printls();
+ DI_Point *midEN2 = quadMidNode(U[0], pt(s3), pt(s1), e, RPNi); // intersection between ls and the bissector between the cutting points
+ cuttingPoints.push_back(new DI_CuttingPoint(midEN2));
+ t1->addQuadEdge (2, midEN2, e, RPNi);
+ t2->addQuadEdge (1, midEN2, e, RPNi);
+ ln->addQuadEdge (0, midEN2, e, RPNi);
+ delete midEN2;
}
subTriangles.push_back(t1);
subTriangles.push_back(t2);
surfLines.push_back(ln);
+ delete U[0];
break;
}
case 2: {// 2 edges cut => 1 triangle + 1 quad => split into 3 triangles
@@ -1514,26 +1526,27 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevel
int s2 = (s1 + 1) % 3;
int s3 = (s1 + 2) % 3;
if(s1 == 0) { // (0,2) => swap U[0] and U[1]
- swap(U[0], U[1]);
+ swap(&U[0], &U[1]);
}
bool useQuads = false; // if true, split the triangle into 1 triangle and 1 quad
- DI_Triangle t1 = DI_Triangle(U[0], pt(s1), U[1], lsTag());
- DI_Line ln = DI_Line(U[0], U[1], LStag);
- DI_Quad q1; DI_Triangle t2, t3;
- if(useQuads) q1 = DI_Quad(U[0], U[1], pt(s2), pt(s3), lsTag());
- else bestQuality(U[0], U[1], pt(s2), pt(s3), t2, t3);
- t2.setLsTag(lsTag()); t3.setLsTag(lsTag());
+ DI_Triangle *t1 = new DI_Triangle(U[0], pt(s1), U[1], lsTag());
+ DI_Line *ln = new DI_Line(U[0], U[1], LStag);
+ DI_Quad *q1; DI_Triangle *t2, *t3;
+ if(useQuads) q1 = new DI_Quad(U[0], U[1], pt(s2), pt(s3), lsTag());
+ else bestQuality(U[0], U[1], pt(s2), pt(s3), &t2, &t3);
+ t2->setLsTag(lsTag()); t3->setLsTag(lsTag());
if(quadT){
- DI_Point midEN2 = quadMidNode(U[0], U[1], pt(s2), e, RPNi); // intersection between ls and the bissector between the cutting points
- cuttingPoints.push_back(DI_CuttingPoint(midEN2));
+ DI_Point *midEN2 = quadMidNode(U[0], U[1], pt(s2), e, RPNi); // intersection between ls and the bissector between the cutting points
+ cuttingPoints.push_back(new DI_CuttingPoint(midEN2));
//midEN2.printls();
- t1.addQuadEdge(2, &midEN2, e, RPNi); //t1.printls();
- ln.addQuadEdge(0, &midEN2, e, RPNi); //ln.printls();
- if(useQuads) q1.addQuadEdge(0, &midEN2, e, RPNi);
- else t2.addQuadEdge(0, &midEN2, e, RPNi); //t2.printls();
+ t1->addQuadEdge(2, midEN2, e, RPNi); //t1.printls();
+ ln->addQuadEdge(0, midEN2, e, RPNi); //ln.printls();
+ if(useQuads) q1->addQuadEdge(0, midEN2, e, RPNi);
+ else t2->addQuadEdge(0, midEN2, e, RPNi); //t2.printls();
+ delete midEN2;
}
subTriangles.push_back(t1);
surfLines.push_back(ln);
@@ -1542,49 +1555,50 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevel
subTriangles.push_back(t2);
subTriangles.push_back(t3);
}
+ delete U[0]; delete U[1];
break;
}
default:
printf("Error : %d edge(s) cut in the triangle (ls : %g %g %g)\n",
- COUNT, pt(0).ls(), pt(1).ls(), pt(2).ls());
+ COUNT, ls(0), ls(1), ls(2));
}
}
// Split a reference Quadrangle into 2 triangles
-void DI_Quad::splitIntoTriangles(std::vector<DI_Triangle> &triangles) const
+void DI_Quad::splitIntoTriangles(std::vector<DI_Triangle *> &triangles) const
{
- triangles.push_back(DI_Triangle(pt(0), pt(1), pt(3), lsTag())); // 013
- triangles.push_back(DI_Triangle(pt(1), pt(2), pt(3), lsTag())); // 123
+ triangles.push_back(new DI_Triangle(pt(0), pt(1), pt(3), lsTag())); // 013
+ triangles.push_back(new DI_Triangle(pt(1), pt(2), pt(3), lsTag())); // 123
}
// Split a reference tetrahedron cut by a level set (defined in a hex) into sub tetrahedra
// and create triangles on the surface of the level set
void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Tetra> &subTetras, std::vector<DI_Triangle> &surfTriangles,
- std::vector<DI_CuttingPoint> &cuttingPoints, std::vector<DI_QualError> &QError) const
+ std::vector<DI_Tetra *> &subTetras, std::vector<DI_Triangle *> &surfTriangles,
+ std::vector<DI_CuttingPoint *> &cuttingPoints, std::vector<DI_QualError *> &QError) const
{
bool quadT = (e->getPolynomialOrder() == 2) ? true : false; // use of quadratic surf triangles and quadratic sub tetrahedra
int tag = RPNi.back()->getTag();
int nbZe = 0, ze[3];
for(int i = 0; i < 4; i++)
- if(pt(i).isOnBorder())
+ if(pt(i)->isOnBorder())
ze[nbZe++] = i;
for(int i = 0; i < nbZe; i++)
- cuttingPoints.push_back(DI_CuttingPoint(pt(ze[i]))); // !! we add these points several times => remove later
+ cuttingPoints.push_back(new DI_CuttingPoint(pt(ze[i]))); // !! we add these points several times => remove later
// case 0 : the tetrahedron is not cut by the level set
if (!isCrossed(pt(0), pt(1)) && !isCrossed(pt(0), pt(2)) && !isCrossed(pt(1), pt(2)) &&
!isCrossed(pt(0), pt(3)) && !isCrossed(pt(1), pt(3)) && !isCrossed(pt(2), pt(3))) {
- subTetras.push_back(*this);
+ subTetras.push_back(new DI_Tetra(*this));
if(nbZe == 3)
- surfTriangles.push_back(DI_Triangle(pt(ze[0]), pt(ze[1]), pt(ze[2]), tag));
+ surfTriangles.push_back(new DI_Triangle(pt(ze[0]), pt(ze[1]), pt(ze[2]), tag));
// !! we add these triangles twice => remove the second later
return;
}
- DI_Point U[4]; // cutting points
- int COUNT = 0; // number of edges cut
+ DI_Point *U[4]; // cutting points
+ int COUNT = 0; // number of edges cut
int IND[4]; // ids of edges cut
// edges nodes and opposite edges nodes
@@ -1594,13 +1608,12 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
for(int i = 0; i < 6; i++){
int n1 = tab[i][0], n2 = tab[i][1];
if (isCrossed(pt(n1), pt(n2))) {
- //printf("edge %d%d cut vls:%g/%g (%g,%g,%g) (%g,%g,%g)\n",n1,n2,vls[n1],vls[n2],x(n1),y(n1),z(n1),x(n2),y(n2),z(n2));
U[COUNT] = Newton(pt(n1), pt(n2), e, RPNi);
IND[COUNT++] = i;
}
}
for(int i = 0; i < COUNT; i++)
- cuttingPoints.push_back(DI_CuttingPoint(U[i]));
+ cuttingPoints.push_back(new DI_CuttingPoint(U[i]));
// Adjust the indices of the nodes in order to have the same pattern for each case with the same number of edges cut;
// compute the position of the middle nodes on the quadratic edges of the sub tetras;
@@ -1614,31 +1627,33 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
int s3 = tab[i1][2];
int s4 = tab[i1][3];
- DI_Tetra t1 = DI_Tetra(U[0], pt(s4), pt(s3), pt(s1));
- DI_Tetra t2 = DI_Tetra(pt(s4), U[0], pt(s3), pt(s2));
- DI_Triangle tr = DI_Triangle(pt(s4), U[0], pt(s3), tag);
+ DI_Tetra *t1 = new DI_Tetra(U[0], pt(s4), pt(s3), pt(s1));
+ DI_Tetra *t2 = new DI_Tetra(pt(s4), U[0], pt(s3), pt(s2));
+ DI_Triangle *tr = new DI_Triangle(pt(s4), U[0], pt(s3), tag);
if(quadT){
- DI_Point midEN2[2]; // intersection between ls and the bissector between the cutting points
+ DI_Point *midEN2[2]; // intersection between ls and the bissector between the cutting points
midEN2[0] = quadMidNode(U[0], pt(s4), pt(s1), e, RPNi);
midEN2[1] = quadMidNode(U[0], pt(s3), pt(s1), e, RPNi);
- tr.addQuadEdge (0, &midEN2[0], e, RPNi);
- t1.addQuadEdge (0, &midEN2[0], e, RPNi);
- t2.addQuadEdge (0, &midEN2[0], e, RPNi);
- tr.addQuadEdge (1, &midEN2[1], e, RPNi);
- t1.addQuadEdge (1, &midEN2[1], e, RPNi);
- t2.addQuadEdge (3, &midEN2[1], e, RPNi);
+ tr->addQuadEdge (0, midEN2[0], e, RPNi);
+ t1->addQuadEdge (0, midEN2[0], e, RPNi);
+ t2->addQuadEdge (0, midEN2[0], e, RPNi);
+ tr->addQuadEdge (1, midEN2[1], e, RPNi);
+ t1->addQuadEdge (1, midEN2[1], e, RPNi);
+ t2->addQuadEdge (3, midEN2[1], e, RPNi);
+ delete midEN2[0]; delete midEN2[1];
}
subTetras.push_back(t1);
subTetras.push_back(t2);
surfTriangles.push_back(tr);
+ delete U[0];
break;
}
case 2 : {// 2 edges cut => 1 tetra + 1 pyramid => split into 3 tetras
int i1 = IND[0], i2 = IND[1];
if((i1 == 0 && (i2 == 2 ||i2 == 3)) || (i1 == 1 && i2 == 4) || ((i1 == 2 || i1 == 3) && i2 == 5)) {
- swap(U[0], U[1]); swap(i1, i2);
+ swap(&U[0], &U[1]); swap(i1, i2);
}
//printf("case 2 : ind : %d,%d\n",i1,i2);
@@ -1647,125 +1662,132 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
int s3 = commonV(tab[i1][2], tab[i1][3], tab[i2][0], tab[i2][1]);
int s4 = commonV(tab[i1][0], tab[i1][1], tab[i2][2], tab[i2][3]);
- DI_Tetra t1 = DI_Tetra(U[0], U[1], pt(s2), pt(s1));
- DI_Triangle tr = DI_Triangle(U[0], U[1], pt(s1), tag);
- DI_Tetra t2, t3;
- int qual = bestQuality(U[1], U[0], pt(s4), pt(s3), pt(s1), t2, t3);
+ DI_Tetra *t1 = new DI_Tetra(U[0], U[1], pt(s2), pt(s1));
+ DI_Triangle *tr = new DI_Triangle(U[0], U[1], pt(s1), tag);
+ DI_Tetra *t2, *t3;
+ int qual = bestQuality(U[1], U[0], pt(s4), pt(s3), pt(s1), &t2, &t3);
if(quadT){
- DI_Point midEN2[3]; // intersection between ls and the bissector between the cutting points
+ DI_Point *midEN2[3]; // intersection between ls and the bissector between the cutting points
midEN2[0] = quadMidNode(U[1], pt(s1), pt(s2), e, RPNi);
midEN2[1] = quadMidNode(U[0], pt(s1), pt(s2), e, RPNi);
midEN2[2] = quadMidNode(U[1], U[0], pt(s4), e, RPNi);
- tr.addQuadEdge(1, &midEN2[0], e, RPNi);
- t1.addQuadEdge(5, &midEN2[0], e, RPNi);
- t2.addQuadEdge(2, &midEN2[0], e, RPNi);
- t3.addQuadEdge(2, &midEN2[0], e, RPNi);
-
- tr.addQuadEdge(2, &midEN2[1], e, RPNi);
- t1.addQuadEdge(2, &midEN2[1], e, RPNi);
- t2.addQuadEdge(5, &midEN2[1], e, RPNi);
-
- tr.addQuadEdge(0, &midEN2[2], e, RPNi);
- if(t2.addQuadEdge(0, &midEN2[2], e, RPNi)) {
- t1.addQuadEdge(0, &midEN2[2], e, RPNi);
+ tr->addQuadEdge(1, midEN2[0], e, RPNi);
+ t1->addQuadEdge(5, midEN2[0], e, RPNi);
+ t2->addQuadEdge(2, midEN2[0], e, RPNi);
+ t3->addQuadEdge(2, midEN2[0], e, RPNi);
+
+ tr->addQuadEdge(2, midEN2[1], e, RPNi);
+ t1->addQuadEdge(2, midEN2[1], e, RPNi);
+ t2->addQuadEdge(5, midEN2[1], e, RPNi);
+
+ tr->addQuadEdge(0, midEN2[2], e, RPNi);
+ if(t2->addQuadEdge(0, midEN2[2], e, RPNi)) {
+ t1->addQuadEdge(0, midEN2[2], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s3), pt(s4), e, RPNi);
- DI_Point quad = middle(mid, midEN2[2], e, RPNi);
+ DI_Point* mid = middle(pt(s3), pt(s4), e, RPNi);
+ DI_Point* quad = middle(mid, midEN2[2], e, RPNi);
if(qual == 1){
- t2.addQuadEdge(1, &quad, e, RPNi);
- t3.addQuadEdge(0, &quad, e, RPNi);
+ t2->addQuadEdge(1, quad, e, RPNi);
+ t3->addQuadEdge(0, quad, e, RPNi);
}
else {
- t2.addQuadEdge(3, &quad, e, RPNi);
- t3.addQuadEdge(1, &quad, e, RPNi);
+ t2->addQuadEdge(3, quad, e, RPNi);
+ t3->addQuadEdge(1, quad, e, RPNi);
}
- if(t2.addQuadEdge(0, &midEN2[2], e, RPNi)) {
- t1.addQuadEdge(0, &midEN2[2], e, RPNi);
+ if(t2->addQuadEdge(0, midEN2[2], e, RPNi)) {
+ t1->addQuadEdge(0, midEN2[2], e, RPNi);
}
else printf("unable to add quadratic edge U0U1 to the subtetrahedra cas 2.\n");
+ delete mid; delete quad;
}
+ delete midEN2[0]; delete midEN2[1]; delete midEN2[2];
}
subTetras.push_back(t1);
subTetras.push_back(t2);
subTetras.push_back(t3);
surfTriangles.push_back(tr);
+ delete U[0]; delete U[1];
break;
}
case 3 : {// 3 edges cut => 1 tetra + 1 prism => split into 4 tetras
int i1 = IND[0], i2 = IND[1], i3 = IND[2];
if(i1 == 0 && i2 == 3) {
- swap(U[1], U[2]); swap(i2, i3);
+ swap(&U[1], &U[2]); swap(i2, i3);
}
//printf("case 3 : ind : %d,%d,%d\n",i1,i2,i3);
int s1 = commonV(tab[i1][0], tab[i1][1], tab[i2][0], tab[i2][1]);
int s2 = commonV(tab[i2][2], tab[i2][3], tab[i3][2], tab[i3][3]);
int s3 = commonV(tab[i1][2], tab[i1][3], tab[i3][2], tab[i3][3]);
- int s4 = commonV(tab[i1][2], tab[i1][3], tab[i2][2], tab[i2][3]); //printf("s: %d %d %d %d\n",s1,s2,s3,s4);
+ int s4 = commonV(tab[i1][2], tab[i1][3], tab[i2][2], tab[i2][3]);
- DI_Tetra t1 = DI_Tetra(pt(s1), U[0], U[1], U[2]);
- DI_Triangle tr = DI_Triangle(U[0], U[1], U[2], tag); //t1.print(); tr.print();
- DI_Tetra t2, t3, t4;
- bestQuality(U[0], U[1], U[2], pt(s2), pt(s3), pt(s4), t2, t3, t4, QError); //t2.print(); t3.print(); t4.print();
+ DI_Tetra *t1 = new DI_Tetra(pt(s1), U[0], U[1], U[2]);
+ DI_Triangle *tr = new DI_Triangle(U[0], U[1], U[2], tag);
+ DI_Tetra *t2, *t3, *t4;
+ bestQuality(U[0], U[1], U[2], pt(s2), pt(s3), pt(s4), &t2, &t3, &t4, QError);
/*t2 = DI_Tetra(U[0],V[0],W[0], x(s2),y(s2),z(s2), U[1],V[1],W[1], U[2],V[2],W[2]);
t3 = DI_Tetra(U[2],V[2],W[2], U[1],V[1],W[1], x(s3),y(s3),z(s3), x(s2),y(s2),z(s2));
t4 = DI_Tetra(U[2],V[2],W[2], x(s3),y(s3),z(s3), x(s4),y(s4),z(s4), x(s2),y(s2),z(s2));*/
if(quadT){
- DI_Point midEN2[3]; // intersection between ls and the bissector between the cutting points
+ DI_Point *midEN2[3]; // intersection between ls and the bissector between the cutting points
midEN2[0] = quadMidNode(U[2], U[1], pt(s3), e, RPNi);
midEN2[1] = quadMidNode(U[1], U[0], pt(s2), e, RPNi);
midEN2[2] = quadMidNode(U[2], U[0], pt(s2), e, RPNi);
- tr.addQuadEdge(1, &midEN2[0], e, RPNi);
- if(t3.addQuadEdge(0, &midEN2[0], e, RPNi)) {
- t1.addQuadEdge(5, &midEN2[0], e, RPNi);
- t2.addQuadEdge(4, &midEN2[0], e, RPNi);
+ tr->addQuadEdge(1, midEN2[0], e, RPNi);
+ if(t3->addQuadEdge(0, midEN2[0], e, RPNi)) {
+ t1->addQuadEdge(5, midEN2[0], e, RPNi);
+ t2->addQuadEdge(4, midEN2[0], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s3), pt(s4), e, RPNi);
- DI_Point quad = middle(mid, midEN2[0], e, RPNi);
- t3.addQuadEdge(1, &quad, e, RPNi);
- t4.addQuadEdge(0, &quad, e, RPNi);
- if(t3.addQuadEdge(0, &midEN2[0], e, RPNi)) {
- t1.addQuadEdge(5, &midEN2[0], e, RPNi);
- t2.addQuadEdge(4, &midEN2[0], e, RPNi);
+ DI_Point *mid = middle(pt(s3), pt(s4), e, RPNi);
+ DI_Point *quad = middle(mid, midEN2[0], e, RPNi);
+ t3->addQuadEdge(1, quad, e, RPNi);
+ t4->addQuadEdge(0, quad, e, RPNi);
+ if(t3->addQuadEdge(0, midEN2[0], e, RPNi)) {
+ t1->addQuadEdge(5, midEN2[0], e, RPNi);
+ t2->addQuadEdge(4, midEN2[0], e, RPNi);
}
else printf("unable to add quadratic edge U1U2 to the subtetrahedra cas 3.\n");
+ delete mid; delete quad;
}
- tr.addQuadEdge(0, &midEN2[1], e, RPNi);
- if(t2.addQuadEdge(1, &midEN2[1], e, RPNi)) {
- t1.addQuadEdge(0, &midEN2[1], e, RPNi);
+ tr->addQuadEdge(0, midEN2[1], e, RPNi);
+ if(t2->addQuadEdge(1, midEN2[1], e, RPNi)) {
+ t1->addQuadEdge(0, midEN2[1], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s2), pt(s3), e, RPNi);
- DI_Point quad = middle(mid, midEN2[1], e, RPNi);
- t2.addQuadEdge(3, &quad, e, RPNi);
- t3.addQuadEdge(5, &quad, e, RPNi);
- if(t2.addQuadEdge(1, &midEN2[1], e, RPNi)) {
- t1.addQuadEdge(0, &midEN2[1], e, RPNi);
+ DI_Point *mid = middle(pt(s2), pt(s3), e, RPNi);
+ DI_Point *quad = middle(mid, midEN2[1], e, RPNi);
+ t2->addQuadEdge(3, quad, e, RPNi);
+ t3->addQuadEdge(5, quad, e, RPNi);
+ if(t2->addQuadEdge(1, midEN2[1], e, RPNi)) {
+ t1->addQuadEdge(0, midEN2[1], e, RPNi);
}
else printf("unable to add quadratic edge U1U2 to the subtetrahedra cas 3.\n");
+ delete mid; delete quad;
}
- tr.addQuadEdge(2, &midEN2[2], e, RPNi);
- if(t2.addQuadEdge(2, &midEN2[2], e, RPNi)) {
- t1.addQuadEdge(2, &midEN2[2], e, RPNi);
+ tr->addQuadEdge(2, midEN2[2], e, RPNi);
+ if(t2->addQuadEdge(2, midEN2[2], e, RPNi)) {
+ t1->addQuadEdge(2, midEN2[2], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s2), pt(s4), e, RPNi);
- DI_Point quad = middle(mid, midEN2[2], e, RPNi);
- t2.addQuadEdge(5, &quad, e, RPNi);
- t3.addQuadEdge(2, &quad, e, RPNi);
- t4.addQuadEdge(2, &quad, e, RPNi);
- if(t2.addQuadEdge(2, &midEN2[2], e, RPNi)) {
- t1.addQuadEdge(2, &midEN2[2], e, RPNi);
+ DI_Point *mid = middle(pt(s2), pt(s4), e, RPNi);
+ DI_Point *quad = middle(mid, midEN2[2], e, RPNi);
+ t2->addQuadEdge(5, quad, e, RPNi);
+ t3->addQuadEdge(2, quad, e, RPNi);
+ t4->addQuadEdge(2, quad, e, RPNi);
+ if(t2->addQuadEdge(2, midEN2[2], e, RPNi)) {
+ t1->addQuadEdge(2, midEN2[2], e, RPNi);
}
else printf("unable to add quadratic edge U0U2 to the subtetrahedra cas 3.\n");
+ delete mid; delete quad;
}
+ delete midEN2[0]; delete midEN2[1]; delete midEN2[2];
}
subTetras.push_back(t1);
@@ -1773,31 +1795,30 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
subTetras.push_back(t3);
subTetras.push_back(t4);
surfTriangles.push_back(tr);
+ delete U[0]; delete U[1]; delete U[2];
break;
}
case 4 : {// 4 edges cut => 2 prisms => split into 6 tetras
int i1 = IND[0], i2 = IND[1], i3 = IND[2], i4 = IND[3];
if(i1 == 0 && i2 == 2) {
- swap(U[0], U[1]); swap(i1, i2);
+ swap(&U[0], &U[1]); swap(i1, i2);
}
else if(i1 == 1 && i2 == 2) {
- swap(U[2], U[3]); swap(i3, i4);
+ swap(&U[2], &U[3]); swap(i3, i4);
}
- //printf("case 4 : ind : %d,%d,%d,%d\n",i1,i2,i3,i4);
int s1 = commonV(tab[i1][0], tab[i1][1], tab[i2][0], tab[i2][1]);
int s2 = commonV(tab[i1][0], tab[i1][1], tab[i2][2], tab[i2][3]);
int s3 = commonV(tab[i1][2], tab[i1][3], tab[i2][0], tab[i2][1]);
- int s4 = commonV(tab[i1][2], tab[i1][3], tab[i2][2], tab[i2][3]); //printf("s: %d %d %d %d\n",s1,s2,s3,s4);
-
- DI_Tetra t1, t2, t3, t4, t5, t6;
- DI_Triangle tr1, tr2;
- bestQuality(U[0], U[1], U[2], U[3], tr1, tr2);
- tr1.setLsTag(tag); tr2.setLsTag(tag);
- bestQuality(pt(s1), U[0], U[1], pt(s4), U[3], U[2], t1, t2, t3, QError);
- bestQuality(pt(s2), U[0], U[3], pt(s3), U[1], U[2], t4, t5, t6, QError);
- //tr1.print(); tr2.print(); t1.print(); t2.print(); t3.print(); t4.print(); t5.print(); t6.print();
+ int s4 = commonV(tab[i1][2], tab[i1][3], tab[i2][2], tab[i2][3]);
+
+ DI_Tetra *t1, *t2, *t3, *t4, *t5, *t6;
+ DI_Triangle *tr1, *tr2;
+ bestQuality(U[0], U[1], U[2], U[3], &tr1, &tr2);
+ tr1->setLsTag(tag); tr2->setLsTag(tag);
+ bestQuality(pt(s1), U[0], U[1], pt(s4), U[3], U[2], &t1, &t2, &t3, QError);
+ bestQuality(pt(s2), U[0], U[3], pt(s3), U[1], U[2], &t4, &t5, &t6, QError);
/*t1 = DI_Tetra(U[0],V[0],W[0], x(s2),y(s2),z(s2), U[1],V[1],W[1], U[3],V[3],W[3]);
t2 = DI_Tetra(U[1],V[1],W[1], U[3],V[3],W[3], x(s2),y(s2),z(s2), U[2],V[2],W[2]);
t3 = DI_Tetra(U[1],V[1],W[1], x(s2),y(s2),z(s2), x(s3),y(s3),z(s3), U[2],V[2],W[2]);
@@ -1808,85 +1829,90 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
tr2 = DI_Triangle(U[1],V[1],W[1], U[2],V[2],W[2], U[3],V[3],W[3]);*/
if(quadT){
- DI_Point midEN2[5]; // intersection between ls and the bissector between the cutting points
+ DI_Point *midEN2[5]; // intersection between ls and the bissector between the cutting points
midEN2[0] = quadMidNode(U[1], U[2], pt(s4), e, RPNi);
midEN2[1] = quadMidNode(U[1], U[0], pt(s2), e, RPNi);
midEN2[2] = quadMidNode(U[0], U[3], pt(s4), e, RPNi);
midEN2[3] = quadMidNode(U[2], U[3], pt(s2), e, RPNi);
midEN2[4] = quadMidNode(U[1], U[3], pt(s2), e, RPNi);
//quadMidNode(U[1], U[3], (pt(s2)+pt(s3))/2, (pt(s4)+pt(s1))/2, e, RPNi);
- tr2.addQuadEdge(0, &midEN2[0], e, RPNi);
- if(t6.addQuadEdge(5, &midEN2[0], e, RPNi)) {
- t2.addQuadEdge(2, &midEN2[0], e, RPNi);
- t3.addQuadEdge(2, &midEN2[0], e, RPNi);
+ tr2->addQuadEdge(0, midEN2[0], e, RPNi);
+ if(t6->addQuadEdge(5, midEN2[0], e, RPNi)) {
+ t2->addQuadEdge(2, midEN2[0], e, RPNi);
+ t3->addQuadEdge(2, midEN2[0], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s1), pt(s4), e, RPNi);
- DI_Point quad = middle(mid, midEN2[0], e, RPNi);
- t4.addQuadEdge(5, &quad, e, RPNi);
- t5.addQuadEdge(3, &quad, e, RPNi);
- t6.addQuadEdge(3, &quad, e, RPNi);
- if(t6.addQuadEdge(5, &midEN2[0], e, RPNi)) {
- t2.addQuadEdge(2, &midEN2[0], e, RPNi);
- t3.addQuadEdge(2, &midEN2[0], e, RPNi);
+ DI_Point *mid = middle(pt(s1), pt(s4), e, RPNi);
+ DI_Point *quad = middle(mid, midEN2[0], e, RPNi);
+ t4->addQuadEdge(5, quad, e, RPNi);
+ t5->addQuadEdge(3, quad, e, RPNi);
+ t6->addQuadEdge(3, quad, e, RPNi);
+ if(t6->addQuadEdge(5, midEN2[0], e, RPNi)) {
+ t2->addQuadEdge(2, midEN2[0], e, RPNi);
+ t3->addQuadEdge(2, midEN2[0], e, RPNi);
}
else printf("unable to add quadratic edge U1U2 to the subtetrahedra cas 4.\n");
+ delete mid; delete quad;
}
- tr1.addQuadEdge(0, &midEN2[1], e, RPNi);
- if(t1.addQuadEdge(1, &midEN2[1], e, RPNi)) {
- t4.addQuadEdge(0, &midEN2[1], e, RPNi);
- t5.addQuadEdge(0, &midEN2[1], e, RPNi);
+ tr1->addQuadEdge(0, midEN2[1], e, RPNi);
+ if(t1->addQuadEdge(1, midEN2[1], e, RPNi)) {
+ t4->addQuadEdge(0, midEN2[1], e, RPNi);
+ t5->addQuadEdge(0, midEN2[1], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s2), pt(s3), e, RPNi);
- DI_Point quad = middle(mid, midEN2[1], e, RPNi);
- t1.addQuadEdge(3, &quad, e, RPNi);
- t2.addQuadEdge(1, &quad, e, RPNi);
- t3.addQuadEdge(0, &quad, e, RPNi);
- if(t1.addQuadEdge(1, &midEN2[1], e, RPNi)) {
- t4.addQuadEdge(0, &midEN2[1], e, RPNi);
- t5.addQuadEdge(0, &midEN2[1], e, RPNi);
+ DI_Point *mid = middle(pt(s2), pt(s3), e, RPNi);
+ DI_Point *quad = middle(mid, midEN2[1], e, RPNi);
+ t1->addQuadEdge(3, quad, e, RPNi);
+ t2->addQuadEdge(1, quad, e, RPNi);
+ t3->addQuadEdge(0, quad, e, RPNi);
+ if(t1->addQuadEdge(1, midEN2[1], e, RPNi)) {
+ t4->addQuadEdge(0, midEN2[1], e, RPNi);
+ t5->addQuadEdge(0, midEN2[1], e, RPNi);
}
else printf("unable to add quadratic edge U0U1 to the subtetrahedra cas 4.\n");
+ delete mid; delete quad;
}
- tr1.addQuadEdge(2, &midEN2[2], e, RPNi);
- if(t5.addQuadEdge(2, &midEN2[2], e, RPNi)) {
- t1.addQuadEdge(2, &midEN2[2], e, RPNi);
+ tr1->addQuadEdge(2, midEN2[2], e, RPNi);
+ if(t5->addQuadEdge(2, midEN2[2], e, RPNi)) {
+ t1->addQuadEdge(2, midEN2[2], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s1), pt(s4), e, RPNi);
- DI_Point quad = middle(mid, midEN2[2], e, RPNi);
- t4.addQuadEdge(2, &quad, e, RPNi);
- t5.addQuadEdge(1, &quad, e, RPNi);
- if(t5.addQuadEdge(2, &midEN2[2], e, RPNi)) {
- t1.addQuadEdge(2, &midEN2[2], e, RPNi);
+ DI_Point *mid = middle(pt(s1), pt(s4), e, RPNi);
+ DI_Point *quad = middle(mid, midEN2[2], e, RPNi);
+ t4->addQuadEdge(2, quad, e, RPNi);
+ t5->addQuadEdge(1, quad, e, RPNi);
+ if(t5->addQuadEdge(2, midEN2[2], e, RPNi)) {
+ t1->addQuadEdge(2, midEN2[2], e, RPNi);
}
else printf("unable to add quadratic edge U0U3 to the subtetrahedra cas 4.\n");
+ delete mid; delete quad;
}
- tr2.addQuadEdge(1, &midEN2[3], e, RPNi);
- if(t2.addQuadEdge(5, &midEN2[3], e, RPNi)) {
- t6.addQuadEdge(2, &midEN2[3], e, RPNi);
+ tr2->addQuadEdge(1, midEN2[3], e, RPNi);
+ if(t2->addQuadEdge(5, midEN2[3], e, RPNi)) {
+ t6->addQuadEdge(2, midEN2[3], e, RPNi);
}
else {
- DI_Point mid = middle(pt(s2), pt(s3), e, RPNi);
- DI_Point quad = middle(mid, midEN2[3], e, RPNi);
- t2.addQuadEdge(4, &quad, e, RPNi);
- t3.addQuadEdge(5, &quad, e, RPNi);
- if(t2.addQuadEdge(5, &midEN2[3], e, RPNi)) {
- t6.addQuadEdge(2, &midEN2[3], e, RPNi);
+ DI_Point *mid = middle(pt(s2), pt(s3), e, RPNi);
+ DI_Point *quad = middle(mid, midEN2[3], e, RPNi);
+ t2->addQuadEdge(4, quad, e, RPNi);
+ t3->addQuadEdge(5, quad, e, RPNi);
+ if(t2->addQuadEdge(5, midEN2[3], e, RPNi)) {
+ t6->addQuadEdge(2, midEN2[3], e, RPNi);
}
else printf("unable to add quadratic edge U2U3 to the subtetrahedra cas 4.\n");
+ delete mid; delete quad;
}
- tr1.addQuadEdge(1, &midEN2[4], e, RPNi);
- tr2.addQuadEdge(2, &midEN2[4], e, RPNi);
- t1.addQuadEdge(4, &midEN2[4], e, RPNi);
- t2.addQuadEdge(0, &midEN2[4], e, RPNi);
- t5.addQuadEdge(5, &midEN2[4], e, RPNi);
- t6.addQuadEdge(0, &midEN2[4], e, RPNi);
+ tr1->addQuadEdge(1, midEN2[4], e, RPNi);
+ tr2->addQuadEdge(2, midEN2[4], e, RPNi);
+ t1->addQuadEdge(4, midEN2[4], e, RPNi);
+ t2->addQuadEdge(0, midEN2[4], e, RPNi);
+ t5->addQuadEdge(5, midEN2[4], e, RPNi);
+ t6->addQuadEdge(0, midEN2[4], e, RPNi);
+ for(int i = 0; i < 5; i++) delete midEN2[i];
}
subTetras.push_back(t1);
@@ -1897,23 +1923,24 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
subTetras.push_back(t6);
surfTriangles.push_back(tr1);
surfTriangles.push_back(tr2);
+ for(int i = 0; i < 4; i++) delete U[i];
break;
}
default:
printf("Error : %d edge(s) cut in the tetrahedron (ls : %g %g %g %g)\n",
- COUNT, pt(0).ls(), pt(1).ls(), pt(2).ls(), pt(3).ls());
+ COUNT, ls(0), ls(1), ls(2), ls(3));
}
}
// Split a reference Hexahedron into 6 tetrahedra
-void DI_Hexa::splitIntoTetras(std::vector<DI_Tetra> &tetras) const
+void DI_Hexa::splitIntoTetras(std::vector<DI_Tetra *> &tetras) const
{
- tetras.push_back(DI_Tetra(pt(0), pt(1), pt(3), pt(4))); // 0134
- tetras.push_back(DI_Tetra(pt(1), pt(4), pt(5), pt(7))); // 1457
- tetras.push_back(DI_Tetra(pt(1), pt(3), pt(4), pt(7))); // 1347
- tetras.push_back(DI_Tetra(pt(2), pt(5), pt(6), pt(7))); // 2567
- tetras.push_back(DI_Tetra(pt(1), pt(2), pt(3), pt(7))); // 1237
- tetras.push_back(DI_Tetra(pt(1), pt(5), pt(2), pt(7))); // 1527
+ tetras.push_back(new DI_Tetra(pt(0), pt(1), pt(3), pt(4))); // 0134
+ tetras.push_back(new DI_Tetra(pt(1), pt(4), pt(5), pt(7))); // 1457
+ tetras.push_back(new DI_Tetra(pt(1), pt(3), pt(4), pt(7))); // 1347
+ tetras.push_back(new DI_Tetra(pt(2), pt(5), pt(6), pt(7))); // 2567
+ tetras.push_back(new DI_Tetra(pt(1), pt(2), pt(3), pt(7))); // 1237
+ tetras.push_back(new DI_Tetra(pt(1), pt(5), pt(2), pt(7))); // 1527
}
// -----------------------------------------------------------------------------
@@ -1921,52 +1948,52 @@ void DI_Hexa::splitIntoTetras(std::vector<DI_Tetra> &tetras) const
// -----------------------------------------------------------------------------
// return the integration points in the reference line
-void DI_Line::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint> &ip) const {
+void DI_Line::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
int N = getNGQLPts(polynomialOrder);
IntPt* intpt (getGQLPts(polynomialOrder));
for (int i = 0; i < N; ++i){
- ip.push_back(DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
+ ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
}
}
// return the integration points in the reference triangle
-void DI_Triangle::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint> &ip) const {
+void DI_Triangle::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
int pO = polynomialOrder;
if(pO == 11 || pO == 16 || pO == 18 || pO == 20) pO++;
if(pO == 15) pO = 17;
int N = getNGQTPts(pO);
IntPt* intpt (getGQTPts(pO));
for (int i = 0; i < N; ++i){
- ip.push_back(DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
+ ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
}
}
// return the integration points in the reference triangle
-void DI_Quad::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint> &ip) const {
+void DI_Quad::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
int N = getNGQQPts(polynomialOrder);
IntPt* intpt (getGQQPts(polynomialOrder));
for (int i = 0; i < N; ++i){
- ip.push_back(DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
+ ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
}
}
// return the integration points in the reference tetra
-void DI_Tetra::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint> &ip) const {
+void DI_Tetra::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
int pO = polynomialOrder;
if(pO == 9) pO++;
int N = getNGQTetPts(pO);
IntPt* intpt (getGQTetPts(pO));
for (int i = 0; i < N; ++i){
- ip.push_back(DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
+ ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
}
}
// return the integration points in the reference Cube
-void DI_Hexa::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint> &ip) const {
+void DI_Hexa::getRefIntegrationPoints ( const int polynomialOrder , std::vector<DI_IntegrationPoint *> &ip) const {
int N = getNGQHPts(polynomialOrder);
IntPt* intpt (getGQHPts(polynomialOrder));
for (int i = 0; i < N; ++i){
- ip.push_back(DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
+ ip.push_back(new DI_IntegrationPoint(intpt[i].pt[0], intpt[i].pt[1], intpt[i].pt[2], intpt[i].weight));
}
}
@@ -1975,679 +2002,709 @@ void DI_Hexa::getRefIntegrationPoints ( const int polynomialOrder , std::vector<
// -----------------------------------------------------------------------------
// cut a line into sublines along the levelset curve
-bool DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_CuttingPoint> &cp, const int polynomialOrder,
- std::vector<DI_Line> &lines, int recurLevel, std::map<int, double> nodeLs[2]) const
+bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+ DI_Point::Container &cp, const int polynomialOrder,
+ std::vector<DI_Line *> &lines, int recurLevel, double **nodeLs) const
{
//printf(" "); print();
- std::vector<const gLevelset *> RPN, RPNi;
- Ls.getRPN(RPN);
+ std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
- DI_Line ll(-1, 0, 0, 1, 0, 0, 2.); //reference line
- //ll.setPolynomialOrder(2);
- std::vector<DI_Line> ll_subLines;
- std::vector<DI_CuttingPoint> ll_cp;
+ DI_Line *ll = new DI_Line(-1, 0, 0, 1, 0, 0, 2.); //reference line
+ ll->setPolynomialOrder(getPolynomialOrder());
+ std::vector<DI_Line *> ll_subLines;
+ std::vector<DI_CuttingPoint *> ll_cp;
- RecurElement re(&ll);
- bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
- pushSubElements(&re, ll_subLines);
+ RecurElement *re = new RecurElement(ll);
+ bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
+ pushSubElements(re, ll_subLines);
+ delete re;
if(signChange){
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- ll.addLs(this, *Lsi, nodeLs);
+ ll->addLs(this, Lsi, l, nodeLs);
int nbSubLn = ll_subLines.size();
int nbCP = ll_cp.size();
- for(int i = 0; i < nbSubLn; i++) ll_subLines[i].addLs(&ll);
- for(int i = 0; i < nbCP; i++) ll_cp[i].addLs(&ll);
+ for(int i = 0; i < nbSubLn; i++) ll_subLines[i]->addLs(ll);
+ for(int i = 0; i < nbCP; i++) ll_cp[i]->addLs(ll);
for(int ln = 0; ln < nbSubLn; ln++){
- DI_Line lnt = ll_subLines[0];
+ DI_Line *lnt = ll_subLines[0];
ll_subLines.erase(ll_subLines.begin());
- lnt.cut(&ll, RPNi, ll_subLines, ll_cp);
+ bool c = lnt->cut(ll, RPNi, ll_subLines, ll_cp);
+ if(c) delete lnt;
}
}
else {
for(int l = 0; l < (int)ll_subLines.size(); l++)
- ll_subLines[l].chooseLs(Lsi);
+ ll_subLines[l]->chooseLs(Lsi);
for(int p = 0; p < (int)ll_cp.size(); p++)
- ll_cp[p].chooseLs(Lsi);
+ ll_cp[p]->chooseLs(Lsi);
}
}
}
else{
- ll_subLines[0].addLs(this, Ls);
+ ll_subLines[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- ll.addLs(this, *Lsi, nodeLs);
+ ll->addLs(this, Lsi, l, nodeLs);
}
}
}
for(int l = 0; l < (int)ll_subLines.size(); l++) {
- ll_subLines[l].computeLsTagDom(&ll, RPN);
- DI_Line ll_subLn = ll_subLines[l];
- mappingEl(&ll_subLn);
- ll_subLn.integrationPoints(polynomialOrder, &ll_subLines[l], &ll, RPN, ip);
- lines.push_back(ll_subLn);
+ ll_subLines[l]->computeLsTagDom(ll, RPN);
+ DI_Line *ll_subLn = new DI_Line(*ll_subLines[l]);
+ mappingEl(ll_subLines[l]);
+ ll_subLines[l]->integrationPoints(polynomialOrder, ll_subLn, ll, RPN, ip);
+ lines.push_back(ll_subLines[l]);
+ delete ll_subLn;
}
for(int p = 0; p < (int)ll_cp.size(); p++) {
- if(ll_cp[p].ls() != 0) continue;
+ if(ll_cp[p]->ls() != 0) {delete ll_cp[p]; continue;}
mappingCP(ll_cp[p]);
- bool isIn = false;
- for(int i = (int)cp.size() - 1; i>= 0; i--)
- if(cp[i].equal(ll_cp[p])) {isIn = true; break;}
- if(!isIn) cp.push_back(ll_cp[p]);
+ std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(ll_cp[p]);
+ if(!isIn.second) delete ll_cp[p];
}
+ delete ll;
return signChange;
}
-// cut a line into sublines along one levelset primitive
-void DI_Line::cut(const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cp) const
+// cut a line into sublines along one levelset primitive and return true if it is cut
+bool DI_Line::cut(const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint *> &cp)
{
// check if the line is cut by the level set
int on = 0, pos = 0, neg = 0, ze[2];
for (int i = 0; i < 2; i++){
- if(pt(i).isOnBorder()) ze[on++] = i;
- else if (pt(i).ls() > 0.) pos++;
+ if(pt(i)->isOnBorder()) ze[on++] = i;
+ else if (ls(i) > 0.) pos++;
else neg++;
}
if(pos && neg)
selfSplit(e, RPNi, subLines, cp);
else {
for(int i = 0; i < on; i++)
- cp.push_back(DI_CuttingPoint(pt(ze[i])));
- subLines.push_back(*this);
+ cp.push_back(new DI_CuttingPoint(pt(ze[i])));
+ subLines.push_back(this);
}
+ return (pos && neg);
}
// cut a triangle into subtriangles along the levelset curve
-bool DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, int recurLevel, std::map<int, double> nodeLs[3]) const
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, int recurLevel, double **nodeLs) const
{
//printf(" ");print();
- std::vector<const gLevelset *> RPN, RPNi;
- Ls.getRPN(RPN);
+ std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
- DI_Triangle tt(0, 0, 0, 1, 0, 0, 0, 1, 0, 0.5); //reference triangle
- //tt.setPolynomialOrder(2);
- std::vector<DI_Quad> tt_subQuads;
- std::vector<DI_Triangle> tt_subTriangles;
- std::vector<DI_Line> tt_surfLines;
- std::vector<DI_CuttingPoint> tt_cp;
+ DI_Triangle *tt = new DI_Triangle(0, 0, 0, 1, 0, 0, 0, 1, 0, 0.5); //reference triangle
+ tt->setPolynomialOrder(getPolynomialOrder());
+ std::vector<DI_Quad *> tt_subQuads;
+ std::vector<DI_Triangle *> tt_subTriangles;
+ std::vector<DI_Line *> tt_surfLines;
+ std::vector<DI_CuttingPoint *> tt_cp;
- RecurElement re(&tt);
- bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
- pushSubElements(&re, tt_subTriangles);
+ RecurElement *re = new RecurElement(tt);
+ bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
+ pushSubElements(re, tt_subTriangles);
+ delete re;
if(signChange){
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l]; //printf("LS(0,0)=%g LS(1,1)=%g\n",(*Lsi)(0,0,0),(*Lsi)(1,1,0));
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- tt.addLs(this, *Lsi, nodeLs);
+ tt->addLs(this, Lsi, l, nodeLs);
int nbSubQ = tt_subQuads.size(), nbSubQ1 = nbSubQ;
int nbSubTr = tt_subTriangles.size(), nbSubTr1 = nbSubTr;
int nbSurfLn = tt_surfLines.size(), nbSurfLn1 = nbSurfLn;
int nbCP = tt_cp.size();
- for(int i = 0; i < nbSubQ; i++) tt_subQuads[i].addLs(&tt);
- for(int i = 0; i < nbSubTr; i++) tt_subTriangles[i].addLs(&tt);
- for(int i = 0; i < nbSurfLn; i++) tt_surfLines[i].addLs(&tt);
- for(int i = 0; i < nbCP; i++) tt_cp[i].addLs(&tt);
+ for(int i = 0; i < nbSubQ; i++) tt_subQuads[i]->addLs(tt);
+ for(int i = 0; i < nbSubTr; i++) tt_subTriangles[i]->addLs(tt);
+ for(int i = 0; i < nbSurfLn; i++) tt_surfLines[i]->addLs(tt);
+ for(int i = 0; i < nbCP; i++) tt_cp[i]->addLs(tt);
int ze[3], cz = 0;
- for (int i = 0; i < 3; i++) if(tt.pt(i).isOnBorder()) ze[cz++] = i;
+ for (int i = 0; i < 3; i++) if(tt->pt(i)->isOnBorder()) ze[cz++] = i;
for(int q = 0; q < nbSubQ; q++) {
nbSubQ1 = tt_subQuads.size();
nbSurfLn1 = tt_surfLines.size();
- DI_Quad qt = tt_subQuads[0];
+ DI_Quad *qt = tt_subQuads[0];
tt_subQuads.erase(tt_subQuads.begin());
- qt.cut(&tt, RPNi, tt_subQuads, tt_subTriangles, tt_surfLines, tt_cp);
+ bool c = qt->cut(tt, RPNi, tt_subQuads, tt_subTriangles, tt_surfLines, tt_cp);
+ if(c) delete qt;
if((int)tt_surfLines.size() - nbSurfLn1 == 1 && (int)tt_subQuads.size() == nbSubQ1){ // 1 line created on boundary of the quad
if (cz == 2) { // 1 line created on boundary of the big triangle => check among surfLines
DI_Line lf (pt(ze[0]), pt(ze[1]));
for(int i = (int)surfLines.size() - 1; i >= 0; i--){
- if (lf.contain(&surfLines[i])){ // line allready created on another surface => pop the new one
- tt_surfLines.pop_back(); break;
+ if (lf.contain(surfLines[i])){ // line allready created on another surface => pop the new one
+ delete tt_surfLines.back(); tt_surfLines.pop_back(); break;
}
}
}
else { // 1 line created inside the big triangle => check among tt_surfLines
- if(isLastLnInV(tt_surfLines))
- tt_surfLines.pop_back();
+ if(isLastLnInV(tt_surfLines)) {
+ delete tt_surfLines.back(); tt_surfLines.pop_back();
+ }
}
}
}
for(int t = 0; t < nbSubTr; t++){
nbSubTr1 = tt_subTriangles.size();
nbSurfLn1 = tt_surfLines.size();
- DI_Triangle trt = tt_subTriangles[0];
+ DI_Triangle *trt = tt_subTriangles[0];
tt_subTriangles.erase(tt_subTriangles.begin());
- trt.cut(&tt, RPNi, tt_subQuads, tt_subTriangles, tt_surfLines, tt_cp);
+ bool c = trt->cut(tt, RPNi, tt_subQuads, tt_subTriangles, tt_surfLines, tt_cp);
+ if(c) delete trt;
if((int)tt_surfLines.size() - nbSurfLn1 == 1 && (int)tt_subTriangles.size() == nbSubTr1) { // 1 line created on boundary of the triangle
if(cz == 2) { // 1 line created on boundary of the big triangle => check among surfLines
DI_Line lf(pt(ze[0]), pt(ze[1]));
for(int i = (int)surfLines.size() - 1; i >= 0; i--){
- if(lf.contain(&surfLines[i])) {
- tt_surfLines.pop_back(); break;
+ if(lf.contain(surfLines[i])) {
+ delete tt_surfLines.back(); tt_surfLines.pop_back(); break;
}
}
}
else{ // 1 line created inside the big triangle => check among tt_surfLines
- if(isLastLnInV(tt_surfLines))
- tt_surfLines.pop_back();
+ if(isLastLnInV(tt_surfLines)) {
+ delete tt_surfLines.back(); tt_surfLines.pop_back();
+ }
}
}
}
for(int ln = 0; ln < nbSurfLn; ln++){
- DI_Line lnt = tt_surfLines[0];
+ DI_Line* lnt = tt_surfLines[0];
tt_surfLines.erase(tt_surfLines.begin());
- lnt.cut(&tt, RPNi, tt_surfLines, tt_cp);
+ bool c = lnt->cut(tt, RPNi, tt_surfLines, tt_cp);
+ if(c) delete lnt;
}
}
else {
for(int q = 0; q < (int)tt_subQuads.size(); q++)
- tt_subQuads[q].chooseLs(Lsi);
+ tt_subQuads[q]->chooseLs(Lsi);
for(int t = 0; t < (int)tt_subTriangles.size(); t++)
- tt_subTriangles[t].chooseLs(Lsi);
+ tt_subTriangles[t]->chooseLs(Lsi);
for(int l = 0; l < (int)tt_surfLines.size(); l++)
- tt_surfLines[l].chooseLs(Lsi);
+ tt_surfLines[l]->chooseLs(Lsi);
for(int p = 0; p < (int)tt_cp.size(); p++)
- tt_cp[p].chooseLs(Lsi);
+ tt_cp[p]->chooseLs(Lsi);
}
}
}
else{
- tt_subTriangles[0].addLs(this, Ls);
+ tt_subTriangles[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- tt.addLs(this, *Lsi, nodeLs);
+ tt->addLs(this, Lsi, l, nodeLs);
}
}
}
for(int q = 0; q < (int)tt_subQuads.size(); q++) {
- tt_subQuads[q].computeLsTagDom(&tt, RPN);
- DI_Quad tt_subQ = tt_subQuads[q];
- mappingEl(&tt_subQ);
- tt_subQ.integrationPoints(polynomialOrderQ, &tt_subQuads[q], &tt, RPN, ip);
- subQuads.push_back(tt_subQ);
+ tt_subQuads[q]->computeLsTagDom(tt, RPN);
+ DI_Quad *tt_subQ = new DI_Quad(*tt_subQuads[q]);
+ mappingEl(tt_subQuads[q]);
+ tt_subQuads[q]->integrationPoints(polynomialOrderQ, tt_subQ, tt, RPN, ip);
+ subQuads.push_back(tt_subQuads[q]);
+ delete tt_subQ;
}
for(int t = 0; t < (int)tt_subTriangles.size(); t++) {
- tt_subTriangles[t].computeLsTagDom(&tt, RPN);
- DI_Triangle tt_subTr = tt_subTriangles[t];
- mappingEl(&tt_subTr);
- tt_subTr.integrationPoints(polynomialOrderTr, &tt_subTriangles[t], &tt, RPN, ip);
- subTriangles.push_back(tt_subTr);
+ tt_subTriangles[t]->computeLsTagDom(tt, RPN);
+ DI_Triangle *tt_subTr = new DI_Triangle(*tt_subTriangles[t]);
+ mappingEl(tt_subTriangles[t]);
+ tt_subTriangles[t]->integrationPoints(polynomialOrderTr, tt_subTr, tt, RPN, ip);
+ subTriangles.push_back(tt_subTriangles[t]);
+ delete tt_subTr;
}
for(int l = 0; l < (int)tt_surfLines.size(); l++) {
- tt_surfLines[l].computeLsTagBound(tt_subQuads, tt_subTriangles);
- if(tt_surfLines[l].lsTag() == -1) continue;
- DI_Line tt_surfLn = tt_surfLines[l];
- mappingEl(&tt_surfLn);
- tt_surfLn.integrationPoints(polynomialOrderL, &tt_surfLines[l], &tt, RPN, ipS);
- surfLines.push_back(tt_surfLn);
+ tt_surfLines[l]->computeLsTagBound(tt_subQuads, tt_subTriangles);
+ if(tt_surfLines[l]->lsTag() == -1) {delete tt_surfLines[l]; continue;}
+ DI_Line *tt_surfLn = new DI_Line(*tt_surfLines[l]);
+ mappingEl(tt_surfLines[l]);
+ tt_surfLines[l]->integrationPoints(polynomialOrderL, tt_surfLn, tt, RPN, ipS);
+ surfLines.push_back(tt_surfLines[l]);
+ delete tt_surfLn;
}
for(int p = 0; p < (int)tt_cp.size(); p++) {
- if(tt_cp[p].ls() != 0) continue;
+ if(tt_cp[p]->ls() != 0) {delete tt_cp[p]; continue;}
mappingCP(tt_cp[p]);
- bool isIn = false;
- for(int i = (int)cp.size() - 1; i>= 0; i--)
- if(cp[i].equal(tt_cp[p])) {isIn = true; break;}
- if(!isIn) cp.push_back(tt_cp[p]);
+ std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(tt_cp[p]);
+ if(!isIn.second) delete tt_cp[p];
}
+ delete tt;
return signChange;
}
// cut a triangle into subtriangles along one levelset primitive
-void DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cp) const
+bool DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cp)
{
// check if the triangle is cut by the level set
int on = 0, pos = 0, neg = 0, ze[3];
for (int i = 0; i < 3; i++){
- if(pt(i).isOnBorder()) ze[on++] = i;
- else if (pt(i).ls() > 0.) pos++;
+ if(pt(i)->isOnBorder()) ze[on++] = i;
+ else if (pt(i)->ls() > 0.) pos++;
else neg++;
}
if(pos && neg)
selfSplit(e, RPNi, subQuads, subTriangles, surfLines, cp);
else {
if(on == 2){ // the level set is zero on an edge of the triangle
- surfLines.push_back(DI_Line(pt(ze[0]), pt(ze[1]), RPNi.back()->getTag()));
+ surfLines.push_back(new DI_Line(pt(ze[0]), pt(ze[1]), RPNi.back()->getTag()));
// add the line twice if the edge belongs to 2 triangles => remove it later!
}
for(int i = 0; i < on; i++)
- cp.push_back(DI_CuttingPoint(pt(ze[i])));
- subTriangles.push_back(*this);
+ cp.push_back(new DI_CuttingPoint(pt(ze[i])));
+ subTriangles.push_back(this);
}
+ return (pos && neg);
}
// cut a quadrangle into subtriangles along the levelset curve
-bool DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, int recurLevel, std::map<int, double> nodeLs[4]) const
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, int recurLevel, double **nodeLs) const
{
- printf(" "); printls();
- std::vector<const gLevelset *> RPN, RPNi;
- Ls.getRPN(RPN);
+ //printf(" "); printls();
+ std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
- DI_Quad qq(-1, -1, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, 4.); // reference quad
- //qq.setPolynomialOrder(2);
- std::vector<DI_Quad> qq_subQuads;
- std::vector<DI_Triangle> qq_subTriangles;
- std::vector<DI_Line> qq_surfLines;
- std::vector<DI_CuttingPoint> qq_cp;
+ DI_Quad *qq = new DI_Quad(-1, -1, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, 4.); // reference quad
+ qq->setPolynomialOrder(getPolynomialOrder());
+ std::vector<DI_Quad *> qq_subQuads;
+ std::vector<DI_Triangle *> qq_subTriangles;
+ std::vector<DI_Line *> qq_surfLines;
+ std::vector<DI_CuttingPoint *> qq_cp;
- RecurElement re(&qq);
- bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
- pushSubElements(&re, qq_subQuads);
+ RecurElement *re = new RecurElement(qq);
+ bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
+ pushSubElements(re, qq_subQuads);
+ delete re;
if(signChange) {
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- qq.addLs(this, *Lsi, nodeLs);
+ qq->addLs(this, Lsi, l, nodeLs);
int nbSubQ = qq_subQuads.size(), nbSubQ1 = nbSubQ;
int nbSubTr = qq_subTriangles.size(), nbSubTr1 = nbSubTr;
int nbSurfLn = qq_surfLines.size(), nbSurfLn1 = nbSurfLn;
int nbCP = qq_cp.size();
- for(int i = 0; i < nbSubQ; i++) qq_subQuads[i].addLs(&qq);
- for(int i = 0; i < nbSubTr; i++) qq_subTriangles[i].addLs(&qq);
- for(int i = 0; i < nbSurfLn; i++) qq_surfLines[i].addLs(&qq);
- for(int i = 0; i < nbCP; i++) qq_cp[i].addLs(&qq);
+ for(int i = 0; i < nbSubQ; i++) qq_subQuads[i]->addLs(qq);
+ for(int i = 0; i < nbSubTr; i++) qq_subTriangles[i]->addLs(qq);
+ for(int i = 0; i < nbSurfLn; i++) qq_surfLines[i]->addLs(qq);
+ for(int i = 0; i < nbCP; i++) qq_cp[i]->addLs(qq);
int pos = 0, neg = 0, ze[2], cz = 0;
for (int i = 0; i < 4; i++){
- if(qq.pt(i).isOnBorder()) ze[cz++] = i;
- else if (qq.ls(i) > 0.) pos++;
+ if(qq->pt(i)->isOnBorder()) ze[cz++] = i;
+ else if (qq->ls(i) > 0.) pos++;
else neg++;
}
for(int q = 0; q < nbSubQ; q++) {
nbSubQ1 = qq_subQuads.size();
nbSurfLn1 = qq_surfLines.size();
- DI_Quad qt = qq_subQuads[0];
+ DI_Quad *qt = qq_subQuads[0];
qq_subQuads.erase(qq_subQuads.begin());
- qt.cut(&qq, RPNi, qq_subQuads, qq_subTriangles, qq_surfLines, qq_cp);
+ bool c = qt->cut(qq, RPNi, qq_subQuads, qq_subTriangles, qq_surfLines, qq_cp);
+ if(c) delete qt;
if((int)qq_surfLines.size() - nbSurfLn1 == 1 && (int)qq_subQuads.size() == nbSubQ1){ // 1 line created on boundary of the quad
if (cz == 2 && !(pos && neg)) { // 1 line created on boundary of the big quad => check among surfLines
DI_Line lf (pt(ze[0]), pt(ze[1]));
for(int i = (int)surfLines.size() - 1; i >= 0; i--){
- if (lf.contain(&surfLines[i])){ // line allready created on another surface => pop the new one
- qq_surfLines.pop_back(); break;
+ if (lf.contain(surfLines[i])){ // line allready created on another surface => pop the new one
+ delete qq_surfLines.back(); qq_surfLines.pop_back(); break;
}
}
}
else { // 1 line created inside the big quad => check among qq_surfLines
- if(isLastLnInV(qq_surfLines))
- qq_surfLines.pop_back();
+ if(isLastLnInV(qq_surfLines)) {
+ delete qq_surfLines.back(); qq_surfLines.pop_back();
+ }
}
}
}
for(int t = 0; t < nbSubTr; t++){
nbSubTr1 = qq_subTriangles.size();
nbSurfLn1 = qq_surfLines.size();
- DI_Triangle trt = qq_subTriangles[0];
+ DI_Triangle *trt = qq_subTriangles[0];
qq_subTriangles.erase(qq_subTriangles.begin());
- trt.cut(&qq, RPNi, qq_subQuads, qq_subTriangles, qq_surfLines, qq_cp);
+ bool c = trt->cut(qq, RPNi, qq_subQuads, qq_subTriangles, qq_surfLines, qq_cp);
+ if(c) delete trt;
if((int)qq_surfLines.size() - nbSurfLn1 == 1 && (int)qq_subTriangles.size() == nbSubTr1) { // 1 line created on boundary of the triangle
if(cz == 2 && !(pos && neg)) { // 1 line created on boundary of the big quad => check among surfLines
DI_Line lf(pt(ze[0]), pt(ze[1]));
for(int i = (int)surfLines.size() - 1; i >= 0; i--){
- if(lf.contain(&surfLines[i])) {
- qq_surfLines.pop_back(); break;
+ if(lf.contain(surfLines[i])) {
+ delete qq_surfLines.back(); qq_surfLines.pop_back(); break;
}
}
}
else{ // 1 line created inside the big quad => check among qq_surfLines
- if(isLastLnInV(qq_surfLines))
- qq_surfLines.pop_back();
+ if(isLastLnInV(qq_surfLines)) {
+ delete qq_surfLines.back(); qq_surfLines.pop_back();
+ }
}
}
}
for(int ln = 0; ln < nbSurfLn; ln++){
- DI_Line lnt = qq_surfLines[0];
+ DI_Line *lnt = qq_surfLines[0];
qq_surfLines.erase(qq_surfLines.begin());
- lnt.cut(&qq, RPNi, qq_surfLines, qq_cp);
+ bool c = lnt->cut(qq, RPNi, qq_surfLines, qq_cp);
+ if(c) delete lnt;
}
}
else {
for(int q = 0; q < (int)qq_subQuads.size(); q++)
- qq_subQuads[q].chooseLs(Lsi);
+ qq_subQuads[q]->chooseLs(Lsi);
for(int t = 0; t < (int)qq_subTriangles.size(); t++)
- qq_subTriangles[t].chooseLs(Lsi);
+ qq_subTriangles[t]->chooseLs(Lsi);
for(int l = 0; l < (int)qq_surfLines.size(); l++)
- qq_surfLines[l].chooseLs(Lsi);
+ qq_surfLines[l]->chooseLs(Lsi);
for(int p = 0; p < (int)qq_cp.size(); p++)
- qq_cp[p].chooseLs(Lsi);
+ qq_cp[p]->chooseLs(Lsi);
}
}
}
else {
- qq_subQuads[0].addLs(this, Ls);
+ qq_subQuads[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- qq.addLs(this, *Lsi, nodeLs);
+ qq->addLs(this, Lsi, l, nodeLs);
}
}
}
for(int q = 0; q < (int)qq_subQuads.size(); q++) {
- qq_subQuads[q].computeLsTagDom(&qq, RPN);
- DI_Quad qq_subQ = qq_subQuads[q];
- mappingEl(&qq_subQ);
- qq_subQ.integrationPoints(polynomialOrderQ, &qq_subQuads[q], &qq, RPN, ip);
- subQuads.push_back(qq_subQ);
+ qq_subQuads[q]->computeLsTagDom(qq, RPN);
+ DI_Quad *qq_subQ = new DI_Quad(*qq_subQuads[q]);
+ mappingEl(qq_subQuads[q]);
+ qq_subQuads[q]->integrationPoints(polynomialOrderQ, qq_subQ, qq, RPN, ip);
+ subQuads.push_back(qq_subQuads[q]);
+ delete qq_subQ;
}
for(int t = 0; t < (int)qq_subTriangles.size(); t++) {
- qq_subTriangles[t].computeLsTagDom(&qq, RPN);
- DI_Triangle qq_subTr = qq_subTriangles[t];
- mappingEl(&qq_subTr); //qq_subTr.printls();
- qq_subTr.integrationPoints(polynomialOrderTr, &qq_subTriangles[t], &qq, RPN, ip);
- subTriangles.push_back(qq_subTr);
+ qq_subTriangles[t]->computeLsTagDom(qq, RPN);
+ DI_Triangle *qq_subTr = new DI_Triangle(*qq_subTriangles[t]);
+ mappingEl(qq_subTriangles[t]);
+ qq_subTriangles[t]->integrationPoints(polynomialOrderTr, qq_subTr, qq, RPN, ip);
+ subTriangles.push_back(qq_subTriangles[t]);
+ delete qq_subTr;
}
for(int l = 0; l < (int)qq_surfLines.size(); l++) {
- qq_surfLines[l].computeLsTagBound(qq_subQuads, qq_subTriangles);
- if(qq_surfLines[l].lsTag() == -1) continue; //FIXME
- DI_Line qq_surfLn = qq_surfLines[l];
- mappingEl(&qq_surfLn);
- qq_surfLn.integrationPoints(polynomialOrderL, &qq_surfLines[l], &qq, RPN, ipS);
- surfLines.push_back(qq_surfLn);
+ qq_surfLines[l]->computeLsTagBound(qq_subQuads, qq_subTriangles);
+ if(qq_surfLines[l]->lsTag() == -1) {delete qq_surfLines[l]; continue;} //FIXME
+ DI_Line *qq_surfLn = new DI_Line(*qq_surfLines[l]);
+ mappingEl(qq_surfLines[l]);
+ qq_surfLines[l]->integrationPoints(polynomialOrderL, qq_surfLn, qq, RPN, ipS);
+ surfLines.push_back(qq_surfLines[l]);
+ delete qq_surfLn;
}
for(int p = 0; p < (int)qq_cp.size(); p++) {
- if(qq_cp[p].ls() != 0) continue;
+ if(qq_cp[p]->ls() != 0) {delete qq_cp[p]; continue;}
mappingCP(qq_cp[p]);
- bool isIn = false;
- for(int i = (int)cp.size() - 1; i >= 0; i--)
- if(cp[i].equal(qq_cp[p])) {isIn = true; break;}
- if(!isIn) cp.push_back(qq_cp[p]);
+ std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(qq_cp[p]);
+ if(!isIn.second) delete qq_cp[p];
}
+ delete qq;
return signChange;
}
// cut a quadrangle into subtriangles along the levelset curve
-void DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cp) const
+bool DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cp)
{
// check if the quadrangle is cut by the level set
int on = 0, pos = 0, neg = 0, ze[4];
for (int i = 0; i < 4; i++){
- if(pt(i).isOnBorder()) ze[on++] = i;
- else if (pt(i).ls() > 0.) pos++;
+ if(pt(i)->isOnBorder()) ze[on++] = i;
+ else if (pt(i)->ls() > 0.) pos++;
else neg++;
}
if (pos && neg) {
- std::vector<DI_Triangle> triangles;
+ std::vector<DI_Triangle *> triangles;
splitIntoTriangles (triangles); // Split the quad into 2 sub triangles
int nl0 = surfLines.size(), nt1, nl1, nt2, nl2;
for (int t = 0; t < (int)triangles.size(); t++) {
nt1 = subTriangles.size(); nl1 = surfLines.size();
- triangles[t].selfSplit(e, RPNi, subQuads, subTriangles, surfLines, cp);
+ triangles[t]->selfSplit(e, RPNi, subQuads, subTriangles, surfLines, cp);
nt2 = subTriangles.size(); nl2 = surfLines.size();
if((nt2 - nt1) == 1 && (nl2 - nl1) == 1){ // only 1 line created on an edge of a triangle => check if it was not yet created
- if(isLastLnInV(surfLines, nl0)) {surfLines.pop_back(); nl2--;}
+ if(isLastLnInV(surfLines, nl0)) {delete surfLines.back(); surfLines.pop_back(); nl2--;}
}
}
}
else {
if(on == 2){ // the level set is zero on an edge of the quad
- surfLines.push_back(DI_Line(pt(ze[0]), pt(ze[1]), RPNi.back()->getTag()));
+ surfLines.push_back(new DI_Line(pt(ze[0]), pt(ze[1]), RPNi.back()->getTag()));
// we add the line twice if the edge belongs to 2 quads => remove it later!
}
for(int i = 0; i < on; i++)
- cp.push_back(DI_CuttingPoint(pt(ze[i])));
- subQuads.push_back(*this);
+ cp.push_back(new DI_CuttingPoint(pt(ze[i])));
+ subQuads.push_back(this);
}
+ return (pos && neg);
}
// cut a tetrahedron into subtetrahedra along the levelset surface
-bool DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderT, const int polynomialOrderQ, const int polynomialOrderTr,
- std::vector<DI_Tetra> &subTetras, std::vector<DI_Quad> &surfQuads,
- std::vector<DI_Triangle> &surfTriangles, int recurLevel, std::map<int, double> nodeLs[4]) const
+ std::vector<DI_Tetra *> &subTetras, std::vector<DI_Quad *> &surfQuads,
+ std::vector<DI_Triangle *> &surfTriangles, int recurLevel, double **nodeLs) const
{
//printf(" ");print();
- std::vector<const gLevelset *> RPN, RPNi;
- Ls.getRPN(RPN);
-
- DI_Tetra tt(0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1./6.); // reference tetra
- //tt.setPolynomialOrder(2);
- std::vector<DI_Hexa> tt_subHexas;
- std::vector<DI_Tetra> tt_subTetras;
- std::vector<DI_Quad> tt_surfQuads;
- std::vector<DI_Triangle> tt_surfTriangles;
- std::vector<DI_Line> tt_surfLines; // not used ...
- std::vector<DI_CuttingPoint> tt_cp;
- std::vector<DI_QualError> QError;
-
- RecurElement re(&tt);
- bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
- pushSubElements(&re, tt_subTetras);
+ std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
+
+ DI_Tetra *tt = new DI_Tetra(0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1./6.); // reference tetra
+ tt->setPolynomialOrder(getPolynomialOrder());
+ std::vector<DI_Hexa *> tt_subHexas;
+ std::vector<DI_Tetra *> tt_subTetras;
+ std::vector<DI_Quad *> tt_surfQuads;
+ std::vector<DI_Triangle *> tt_surfTriangles;
+ std::vector<DI_Line *> tt_surfLines; // not used ...
+ std::vector<DI_CuttingPoint *> tt_cp;
+ std::vector<DI_QualError *> QError;
+
+ RecurElement *re = new RecurElement(tt);
+ bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
+ pushSubElements(re, tt_subTetras);
+ delete re;
if(signChange) {
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- tt.addLs(this, *Lsi, nodeLs);
+ tt->addLs(this, Lsi, l, nodeLs);
int nbSubT = tt_subTetras.size(), nbSubT1 = nbSubT;
int nbSurfQ = tt_surfQuads.size();
int nbSurfTr = tt_surfTriangles.size(), nbSurfTr1 = nbSurfTr;
int nbCP = tt_cp.size();
- for(int i = 0; i < nbSubT; i++) tt_subTetras[i].addLs(&tt);
- for(int i = 0; i < nbSurfQ; i++) tt_surfQuads[i].addLs(&tt);
- for(int i = 0; i < nbSurfTr; i++) tt_surfTriangles[i].addLs(&tt);
- for(int i = 0; i < nbCP; i++) tt_cp[i].addLs(&tt);
+ for(int i = 0; i < nbSubT; i++) tt_subTetras[i]->addLs(tt);
+ for(int i = 0; i < nbSurfQ; i++) tt_surfQuads[i]->addLs(tt);
+ for(int i = 0; i < nbSurfTr; i++) tt_surfTriangles[i]->addLs(tt);
+ for(int i = 0; i < nbCP; i++) tt_cp[i]->addLs(tt);
int ze[3], cz = 0;
- for (int i = 0; i < 4; i++) if(tt.pt(i).isOnBorder()) ze[cz++] = i;
+ for (int i = 0; i < 4; i++) if(tt->pt(i)->isOnBorder()) ze[cz++] = i;
for(int t = 0; t < nbSubT; t++) {
nbSurfTr1 = tt_surfTriangles.size();
nbSubT1 = tt_subTetras.size();
- DI_Tetra tet = tt_subTetras[0];
+ DI_Tetra *tet = tt_subTetras[0];
tt_subTetras.erase(tt_subTetras.begin());
- tet.cut(&tt, RPNi, tt_subTetras, tt_surfQuads, tt_surfTriangles, tt_cp, QError);
+ bool c = tet->cut(tt, RPNi, tt_subTetras, tt_surfQuads, tt_surfTriangles, tt_cp, QError);
+ if(c) delete tet;
if((int)tt_surfTriangles.size() - nbSurfTr1 == 1 && (int)tt_subTetras.size() == nbSubT1) { // 1 triangle created on surface of the subtetra
// check among the tt_surfTriangles
- if(isLastTrInV(tt_surfTriangles))
- tt_surfTriangles.pop_back();
+ if(isLastTrInV(tt_surfTriangles)) {
+ delete tt_surfTriangles.back(); tt_surfTriangles.pop_back();
+ }
else if(cz == 3) { // 1 triangle created on surface of the reference tet => check among surfTriangles
DI_Triangle tf(pt(ze[0]), pt(ze[1]), pt(ze[2]));
for(int i = (int)surfTriangles.size() - 1; i >= 0; i--){
- if(tf.contain(&surfTriangles[i])) {
- tt_surfTriangles.pop_back(); break;
+ if(tf.contain(surfTriangles[i])) {
+ delete tt_surfTriangles.back(); tt_surfTriangles.pop_back(); break;
}
}
}
}
}
for(int q = 0; q < nbSurfQ; q++) {
- DI_Quad qt = tt_surfQuads[0];
+ DI_Quad *qt = tt_surfQuads[0];
tt_surfQuads.erase(tt_surfQuads.begin());
- qt.cut(&tt, RPNi, tt_surfQuads, tt_surfTriangles, tt_surfLines, tt_cp);
+ bool c = qt->cut(tt, RPNi, tt_surfQuads, tt_surfTriangles, tt_surfLines, tt_cp);
+ if(c) delete qt;
}
for(int t = 0; t < nbSurfTr; t++){
- DI_Triangle trt = tt_surfTriangles[0];
+ DI_Triangle *trt = tt_surfTriangles[0];
tt_surfTriangles.erase(tt_surfTriangles.begin());
- trt.cut(&tt, RPNi, tt_surfQuads, tt_surfTriangles, tt_surfLines, tt_cp);
+ bool c = trt->cut(tt, RPNi, tt_surfQuads, tt_surfTriangles, tt_surfLines, tt_cp);
+ if(c) delete trt;
}
}
else {
for(int t = 0; t < (int)tt_subTetras.size(); t++)
- tt_subTetras[t].chooseLs(Lsi);
+ tt_subTetras[t]->chooseLs(Lsi);
for(int q = 0; q < (int)tt_surfQuads.size(); q++)
- tt_surfQuads[q].chooseLs(Lsi);
+ tt_surfQuads[q]->chooseLs(Lsi);
for(int t = 0; t < (int)tt_surfTriangles.size(); t++)
- tt_surfTriangles[t].chooseLs(Lsi);
+ tt_surfTriangles[t]->chooseLs(Lsi);
for(int p = 0; p < (int)tt_cp.size(); p++)
- tt_cp[p].chooseLs(Lsi);
+ tt_cp[p]->chooseLs(Lsi);
}
}
}
else {
- tt_subTetras[0].addLs(this, Ls);
+ tt_subTetras[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- tt.addLs(this, *Lsi, nodeLs);
+ tt->addLs(this, Lsi, l, nodeLs);
}
}
}
- for(int i = 0; i < (int)QError.size(); i++)
- QError[i].print(this);
+ for(int i = 0; i < (int)QError.size(); i++) {
+ QError[i]->print(this);
+ delete QError[i];
+ }
for(int t = 0; t < (int)tt_subTetras.size(); t++) {
- tt_subTetras[t].computeLsTagDom(&tt, RPN);
- DI_Tetra tt_subT = tt_subTetras[t];
- mappingEl(&tt_subT);
- tt_subT.integrationPoints(polynomialOrderT, &tt_subTetras[t], &tt, RPN, ip);
- subTetras.push_back(tt_subT);
+ tt_subTetras[t]->computeLsTagDom(tt, RPN);
+ DI_Tetra *tt_subT = new DI_Tetra(*tt_subTetras[t]);
+ mappingEl(tt_subTetras[t]);
+ tt_subTetras[t]->integrationPoints(polynomialOrderT, tt_subT, tt, RPN, ip);
+ subTetras.push_back(tt_subTetras[t]);
+ delete tt_subT;
}
for(int q = 0; q < (int)tt_surfQuads.size(); q++) {
- tt_surfQuads[q].computeLsTagBound(tt_subHexas, tt_subTetras);
- if(tt_surfQuads[q].lsTag() == -1) continue;
- DI_Quad tt_surfQ = tt_surfQuads[q];
- mappingEl(&tt_surfQ);
- tt_surfQ.integrationPoints(polynomialOrderQ, &tt_surfQuads[q], &tt, RPN, ipS);
- surfQuads.push_back(tt_surfQ);
+ tt_surfQuads[q]->computeLsTagBound(tt_subHexas, tt_subTetras);
+ if(tt_surfQuads[q]->lsTag() == -1) {delete tt_surfQuads[q]; continue;}
+ DI_Quad *tt_surfQ = new DI_Quad(*tt_surfQuads[q]);
+ mappingEl(tt_surfQuads[q]);
+ tt_surfQuads[q]->integrationPoints(polynomialOrderQ, tt_surfQ, tt, RPN, ipS);
+ surfQuads.push_back(tt_surfQuads[q]);
+ delete tt_surfQ;
}
for(int t = 0; t < (int)tt_surfTriangles.size(); t++) {
- tt_surfTriangles[t].computeLsTagBound(tt_subHexas, tt_subTetras);
- if(tt_surfTriangles[t].lsTag() == -1) continue;
- DI_Triangle tt_surfTr = tt_surfTriangles[t];
- mappingEl(&tt_surfTriangles[t]);
- tt_surfTriangles[t].integrationPoints(polynomialOrderTr, &tt_surfTr, &tt, RPN, ipS);
+ tt_surfTriangles[t]->computeLsTagBound(tt_subHexas, tt_subTetras);
+ if(tt_surfTriangles[t]->lsTag() == -1) {delete tt_surfTriangles[t]; continue;}
+ DI_Triangle *tt_surfTr = new DI_Triangle(*tt_surfTriangles[t]);
+ mappingEl(tt_surfTriangles[t]);
+ tt_surfTriangles[t]->integrationPoints(polynomialOrderTr, tt_surfTr, tt, RPN, ipS);
surfTriangles.push_back(tt_surfTriangles[t]);
+ delete tt_surfTr;
}
for(int p = 0; p < (int)tt_cp.size(); p++) {
- if(tt_cp[p].ls() != 0) continue;
+ if(tt_cp[p]->ls() != 0) {delete tt_cp[p]; continue;}
mappingCP(tt_cp[p]);
- bool isIn = false;
- for(int i = (int)cp.size() - 1; i >= 0; i--)
- if(cp[i].equal(tt_cp[p])) {isIn = true; break;}
- if(!isIn) cp.push_back(tt_cp[p]);
+ std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(tt_cp[p]);
+ if(!isIn.second) delete tt_cp[p];
}
+ delete tt;
return signChange;
}
// cut a tetrahedron into subtetrahedra along the levelset surface
-void DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Tetra> &subTetras, std::vector<DI_Quad> &surfQuads,
- std::vector<DI_Triangle> &surfTriangles, std::vector<DI_CuttingPoint> &cp,
- std::vector<DI_QualError> &QError) const
+bool DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Tetra *> &subTetras, std::vector<DI_Quad *> &surfQuads,
+ std::vector<DI_Triangle *> &surfTriangles, std::vector<DI_CuttingPoint *> &cp,
+ std::vector<DI_QualError *> &QError)
{
// check if the tetra is cut by the level set
int on = 0, pos = 0, neg = 0, ze[4];
for (int i = 0; i < 4; i++){
- if(pt(i).isOnBorder()) ze[on++] = i;
- else if (pt(i).ls() > 0.) pos++;
+ if(pt(i)->isOnBorder()) ze[on++] = i;
+ else if (pt(i)->ls() > 0.) pos++;
else neg++;
}
if (pos && neg)
selfSplit(e, RPNi, subTetras, surfTriangles, cp, QError);
else{
if(on == 3){ // the level set is zero on a face of the tetra
- surfTriangles.push_back(DI_Triangle(pt(ze[0]), pt(ze[1]), pt(ze[2]), RPNi.back()->getTag())); // add the triangle twice if the face belongs to 2 tetras => remove it later!
+ surfTriangles.push_back(new DI_Triangle(pt(ze[0]), pt(ze[1]), pt(ze[2]), RPNi.back()->getTag()));
+ // add the triangle twice if the face belongs to 2 tetras => remove it later!
}
for(int i = 0; i < on; i++)
- cp.push_back(DI_CuttingPoint(pt(ze[i])));
- subTetras.push_back(*this); // add the tetra to subTetras if it is not cut
+ cp.push_back(new DI_CuttingPoint(pt(ze[i])));
+ subTetras.push_back(this); // add the tetra to subTetras if it is not cut
}
+ return (pos && neg);
}
// cut a hex into subtetras along the levelset surface
-bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderH, const int polynomialOrderT,
const int polynomialOrderQ, const int polynomialOrderTr,
- std::vector<DI_Hexa> &subHexas, std::vector<DI_Tetra> &subTetras,
- std::vector<DI_Quad> &surfQuads, std::vector<DI_Triangle> &surfTriangles,
- std::vector<DI_Line> &frontLines, int recurLevel, std::map<int, double> nodeLs[8]) const
+ std::vector<DI_Hexa *> &subHexas, std::vector<DI_Tetra *> &subTetras,
+ std::vector<DI_Quad *> &surfQuads, std::vector<DI_Triangle *> &surfTriangles,
+ std::vector<DI_Line *> &frontLines, int recurLevel, double **nodeLs) const
{
- printf(" "); print();
-
- std::vector<const gLevelset *> RPN, RPNi;
- Ls.getRPN(RPN);
+ //printf(" "); print();
+ std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
// reference hexa
- DI_Hexa hh(-1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, 8.);
- std::vector<DI_Hexa> hh_subHexas;
- std::vector<DI_Tetra> hh_subTetras;
- std::vector<DI_Quad> hh_surfQuads;
- std::vector<DI_Triangle> hh_surfTriangles;
- std::vector<DI_Line> hh_frontLines;
- std::vector<DI_CuttingPoint> hh_cp;
- std::vector<DI_QualError> QError;
-
- RecurElement re(&hh);
- bool signChange = re.cut(recurLevel, this, Ls, -1., nodeLs);
- pushSubElements(&re, hh_subHexas);
+ DI_Hexa *hh = new DI_Hexa(-1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, -1, 1, 1, 8.);
+ hh->setPolynomialOrder(getPolynomialOrder());
+ std::vector<DI_Hexa *> hh_subHexas;
+ std::vector<DI_Tetra *> hh_subTetras;
+ std::vector<DI_Quad *> hh_surfQuads;
+ std::vector<DI_Triangle *> hh_surfTriangles;
+ std::vector<DI_Line *> hh_frontLines;
+ std::vector<DI_CuttingPoint *> hh_cp;
+ std::vector<DI_QualError *> QError;
+
+ RecurElement *re = new RecurElement(hh);
+ bool signChange = re->cut(recurLevel, this, RPN, -1., nodeLs);
+ pushSubElements(re, hh_subHexas);
+ delete re;
if(signChange){
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- hh.addLs(this, *Lsi, nodeLs);
+ hh->addLs(this, Lsi, l, nodeLs);
int nbSubH = hh_subHexas.size();
int nbSubT = hh_subTetras.size(), nbSubT1 = nbSubT;
int nbSurfQ = hh_surfQuads.size(), nbSurfQ1 = nbSurfQ;
int nbSurfTr = hh_surfTriangles.size(), nbSurfTr1 = nbSurfTr;
int nbFrontLn = hh_frontLines.size();
int nbCP = hh_cp.size();
- for(int i = 0; i < nbSubH; i++) hh_subHexas[i].addLs(&hh);
- for(int i = 0; i < nbSubT; i++) hh_subTetras[i].addLs(&hh);
- for(int i = 0; i < nbSurfQ; i++) hh_surfQuads[i].addLs(&hh);
- for(int i = 0; i < nbSurfTr; i++) hh_surfTriangles[i].addLs(&hh);
- for(int i = 0; i < nbFrontLn; i++) hh_frontLines[i].addLs(&hh);
- for(int i = 0; i < nbCP; i++) hh_cp[i].addLs(&hh);
+ for(int i = 0; i < nbSubH; i++) hh_subHexas[i]->addLs(hh);
+ for(int i = 0; i < nbSubT; i++) hh_subTetras[i]->addLs(hh);
+ for(int i = 0; i < nbSurfQ; i++) hh_surfQuads[i]->addLs(hh);
+ for(int i = 0; i < nbSurfTr; i++) hh_surfTriangles[i]->addLs(hh);
+ for(int i = 0; i < nbFrontLn; i++) hh_frontLines[i]->addLs(hh);
+ for(int i = 0; i < nbCP; i++) hh_cp[i]->addLs(hh);
int pos = 0, neg = 0, ze[4], cz = 0;
for (int i = 0; i < 8; i++){
- if(hh.pt(i).isOnBorder()) ze[cz++] = i;
- else if (hh.ls(i) > 0.) pos++;
+ if(hh->pt(i)->isOnBorder()) ze[cz++] = i;
+ else if (hh->ls(i) > 0.) pos++;
else neg++;
}
for(int h = 0; h < nbSubH; h++) {
nbSurfQ1 = hh_surfQuads.size();
- DI_Hexa ht = hh_subHexas[0];
+ DI_Hexa *ht = hh_subHexas[0];
hh_subHexas.erase(hh_subHexas.begin());
- ht.cut(&hh, RPNi, hh_subHexas, hh_subTetras, hh_surfQuads, hh_surfTriangles, hh_cp, QError);
+ bool c = ht->cut(hh, RPNi, hh_subHexas, hh_subTetras, hh_surfQuads, hh_surfTriangles, hh_cp, QError);
+ if(c) delete ht;
if((int)hh_surfQuads.size() > nbSurfQ1){ // 1 quad created on surface of the subHexa
- if(isLastQInV(hh_surfQuads))
- hh_surfQuads.pop_back();
+ if(isLastQInV(hh_surfQuads)) {
+ delete hh_surfQuads.back(); hh_surfQuads.pop_back();
+ }
else if (cz == 4) {
if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3])))
swap(ze[2], ze[3]);
DI_Quad qf(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]));
for(int i = (int)surfQuads.size() - 1; i >= 0; i--){
- if (qf.contain(&surfQuads[i])){
- hh_surfQuads.pop_back(); break;
+ if (qf.contain(surfQuads[i])){
+ delete hh_surfQuads.back(); hh_surfQuads.pop_back(); break;
}
}
for(int i = (int)surfTriangles.size() - 1; i >= 0; i--){
- if (qf.contain(&surfTriangles[i])) {
- hh_surfQuads.pop_back(); break;
+ if (qf.contain(surfTriangles[i])) {
+ delete hh_surfQuads.back(); hh_surfQuads.pop_back(); break;
}
}
}
@@ -2656,20 +2713,22 @@ bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
for(int t = 0; t < nbSubT; t++) {
nbSurfTr1 = hh_surfTriangles.size();
nbSubT1 = hh_subTetras.size();
- DI_Tetra tet = hh_subTetras[0];
+ DI_Tetra *tet = hh_subTetras[0];
hh_subTetras.erase(hh_subTetras.begin());
- tet.cut(&hh, RPNi, hh_subTetras, hh_surfQuads, hh_surfTriangles, hh_cp, QError);
+ bool c = tet->cut(hh, RPNi, hh_subTetras, hh_surfQuads, hh_surfTriangles, hh_cp, QError);
+ if(c) delete tet;
if((int)hh_surfTriangles.size() - nbSurfTr1 == 1 && (int)hh_subTetras.size() == nbSubT1) { // 1 triangle created on surface of the subtetra
// check among hh_surfTriangles
- if(isLastTrInV(hh_surfTriangles))
- hh_surfTriangles.pop_back();
+ if(isLastTrInV(hh_surfTriangles)) {
+ delete hh_surfTriangles.back(); hh_surfTriangles.pop_back();
+ }
else if (cz == 4 && !(pos && neg)){ // 1 triangle created on surface of the hex => check among surfTriangles
if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3])))
swap(ze[2], ze[3]);
DI_Quad qt(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]));
for(int i = (int)surfTriangles.size() - 1; i >= 0; i--){
- if(qt.contain(&surfTriangles[i])) {
- hh_surfTriangles.pop_back(); break;
+ if(qt.contain(surfTriangles[i])) {
+ delete hh_surfTriangles.back(); hh_surfTriangles.pop_back(); break;
}
}
}
@@ -2677,7 +2736,7 @@ bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
// check among the quads with quality error
for(int i = 0; i < (int)QError.size(); i++){
if(isInQE(hh_surfTriangles[nbSurfTr1], QError[i])) {
- hh_surfTriangles.pop_back();
+ delete hh_surfTriangles.back(); hh_surfTriangles.pop_back();
break;
}
}
@@ -2685,126 +2744,135 @@ bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
}
}
for(int q = 0; q < nbSurfQ; q++) {
- DI_Quad qt = hh_surfQuads[0];
+ DI_Quad *qt = hh_surfQuads[0];
hh_surfQuads.erase(hh_surfQuads.begin());
- qt.cut(&hh, RPNi, hh_surfQuads, hh_surfTriangles, hh_frontLines, hh_cp);
+ bool c = qt->cut(hh, RPNi, hh_surfQuads, hh_surfTriangles, hh_frontLines, hh_cp);
+ if(c) delete qt;
}
for(int t = 0; t < nbSurfTr; t++){
- DI_Triangle trt = hh_surfTriangles[0];
+ DI_Triangle *trt = hh_surfTriangles[0];
hh_surfTriangles.erase(hh_surfTriangles.begin());
- trt.cut(&hh, RPNi, hh_surfQuads, hh_surfTriangles, hh_frontLines, hh_cp);
+ bool c = trt->cut(hh, RPNi, hh_surfQuads, hh_surfTriangles, hh_frontLines, hh_cp);
+ if(c) delete trt;
}
for(int l = 0; l < nbFrontLn; l++){
- DI_Line lnt = hh_frontLines[0];
+ DI_Line *lnt = hh_frontLines[0];
hh_frontLines.erase(hh_frontLines.begin());
- lnt.cut(&hh, RPNi, hh_frontLines, hh_cp);
+ bool c = lnt->cut(hh, RPNi, hh_frontLines, hh_cp);
+ if(c) delete lnt;
}
}
else {
for(int h = 0; h < (int)hh_subHexas.size(); h++)
- hh_subHexas[h].chooseLs(Lsi);
+ hh_subHexas[h]->chooseLs(Lsi);
for(int t = 0; t < (int)hh_subTetras.size(); t++)
- hh_subTetras[t].chooseLs(Lsi);
+ hh_subTetras[t]->chooseLs(Lsi);
for(int q = 0; q < (int)hh_surfQuads.size(); q++)
- hh_surfQuads[q].chooseLs(Lsi);
+ hh_surfQuads[q]->chooseLs(Lsi);
for(int t = 0; t < (int)hh_surfTriangles.size(); t++)
- hh_surfTriangles[t].chooseLs(Lsi);
+ hh_surfTriangles[t]->chooseLs(Lsi);
for(int l = 0; l < (int)hh_frontLines.size(); l++)
- hh_frontLines[l].chooseLs(Lsi);
+ hh_frontLines[l]->chooseLs(Lsi);
for(int p = 0; p < (int)hh_cp.size(); p++)
- hh_cp[p].chooseLs(Lsi);
+ hh_cp[p]->chooseLs(Lsi);
}
//printf("Lsi=%d sizeLs=%d nbH=%d nbT=%d nbQ=%d nbTr=%d sizeLs(hh_cp)=%d\n",Lsi->type(),(hh_subHexas.size()>0)?hh_subHexas[0].sizeLs():hh_subTetras[0].sizeLs(), hh_subHexas.size(),hh_subTetras.size(),hh_surfQuads.size(),hh_surfTriangles.size(),(hh_cp.size() > 0)?hh_cp[0].sizeLs():0);
}
}
else {
- hh_subHexas[0].addLs(this, Ls);
+ hh_subHexas[0]->addLs(this, RPN.back());
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- hh.addLs(this, *Lsi, nodeLs);
+ hh->addLs(this, Lsi, l, nodeLs);
}
}
}
- for(int i = 0; i < (int)QError.size(); i++)
- QError[i].print(this);
+ for(int i = 0; i < (int)QError.size(); i++) {
+ QError[i]->print(this);
+ delete QError[i];
+ }
for(int h = 0; h < (int)hh_subHexas.size(); h++) {
- hh_subHexas[h].computeLsTagDom(&hh, RPN);
- DI_Hexa hh_subH = hh_subHexas[h];
- mappingEl(&hh_subH);
- hh_subH.integrationPoints(polynomialOrderH, &hh_subHexas[h], &hh, RPN, ip);
- subHexas.push_back(hh_subH);
+ hh_subHexas[h]->computeLsTagDom(hh, RPN);
+ DI_Hexa *hh_subH = new DI_Hexa(*hh_subHexas[h]);
+ mappingEl(hh_subHexas[h]);
+ hh_subHexas[h]->integrationPoints(polynomialOrderH, hh_subH, hh, RPN, ip);
+ subHexas.push_back(hh_subHexas[h]);
+ delete hh_subH;
}
for(int t = 0; t < (int)hh_subTetras.size(); t++) {
- hh_subTetras[t].computeLsTagDom(&hh, RPN);
- DI_Tetra hh_subT = hh_subTetras[t];
- mappingEl(&hh_subT);
- hh_subT.integrationPoints(polynomialOrderT, &hh_subTetras[t], &hh, RPN, ip);
- subTetras.push_back(hh_subT);
+ hh_subTetras[t]->computeLsTagDom(hh, RPN);
+ DI_Tetra *hh_subT = new DI_Tetra(*hh_subTetras[t]);
+ mappingEl(hh_subTetras[t]);
+ hh_subTetras[t]->integrationPoints(polynomialOrderT, hh_subT, hh, RPN, ip);
+ subTetras.push_back(hh_subTetras[t]);
+ delete hh_subT;
}
for(int q = 0; q < (int)hh_surfQuads.size(); q++) {
- hh_surfQuads[q].computeLsTagBound(hh_subHexas, hh_subTetras);
- if(hh_surfQuads[q].lsTag() == -1) continue;
- DI_Quad hh_surfQ = hh_surfQuads[q];
- mappingEl(&hh_surfQ);
- hh_surfQ.integrationPoints(polynomialOrderQ, &hh_surfQuads[q], &hh, RPN, ipS);
- surfQuads.push_back(hh_surfQ);
+ hh_surfQuads[q]->computeLsTagBound(hh_subHexas, hh_subTetras);
+ if(hh_surfQuads[q]->lsTag() == -1) {delete hh_surfQuads[q]; continue;}
+ DI_Quad *hh_surfQ = new DI_Quad(*hh_surfQuads[q]);
+ mappingEl(hh_surfQuads[q]);
+ hh_surfQuads[q]->integrationPoints(polynomialOrderQ, hh_surfQ, hh, RPN, ipS);
+ surfQuads.push_back(hh_surfQuads[q]);
+ delete hh_surfQ;
}
for(int t = 0; t < (int)hh_surfTriangles.size(); t++) {
- hh_surfTriangles[t].computeLsTagBound(hh_subHexas, hh_subTetras);
- if(hh_surfTriangles[t].lsTag() == -1) continue;
- DI_Triangle hh_surfTr = hh_surfTriangles[t];
- mappingEl(&hh_surfTr);
- hh_surfTr.integrationPoints(polynomialOrderTr, &hh_surfTriangles[t], &hh, RPN, ipS);
- surfTriangles.push_back(hh_surfTr);
+ hh_surfTriangles[t]->computeLsTagBound(hh_subHexas, hh_subTetras);
+ if(hh_surfTriangles[t]->lsTag() == -1) {delete hh_surfTriangles[t]; continue;}
+ DI_Triangle *hh_surfTr = new DI_Triangle(*hh_surfTriangles[t]);
+ mappingEl(hh_surfTriangles[t]);
+ hh_surfTriangles[t]->integrationPoints(polynomialOrderTr, hh_surfTr, hh, RPN, ipS);
+ surfTriangles.push_back(hh_surfTriangles[t]);
+ delete hh_surfTr;
}
for(int l = 0; l < (int)hh_frontLines.size(); l++) {
- hh_frontLines[l].computeLsTagBound(hh_surfQuads, hh_surfTriangles);
- if(hh_frontLines[l].lsTag() == -1) continue;
- DI_Line hh_frontLn = hh_frontLines[l];
- mappingEl(&hh_frontLn); //tt_surfLn.printls();
- //hh_frontLn.integrationPoints(polynomialOrderL, &tt_surfLines[l], &tt, RPN, ipS);
- frontLines.push_back(hh_frontLn);
+ hh_frontLines[l]->computeLsTagBound(hh_surfQuads, hh_surfTriangles);
+ if(hh_frontLines[l]->lsTag() == -1) {delete hh_frontLines[l]; continue;}
+ DI_Line *hh_frontLn = new DI_Line(*hh_frontLines[l]);
+ mappingEl(hh_frontLines[l]);
+ //hh_frontLines[l]->integrationPoints(polynomialOrderL, tt_surfLn, tt, RPN, ipS);
+ frontLines.push_back(hh_frontLines[l]);
+ delete hh_frontLn;
}
for(int p = 0; p < (int)hh_cp.size(); p++) {
- if(hh_cp[p].ls() != 0) continue; // returns only the cutting points with ls==0
+ if(hh_cp[p]->ls() != 0) {delete hh_cp[p]; continue;} // returns only the cutting points with ls==0
mappingCP(hh_cp[p]);
- bool isIn = false;
- for(int i = (int)cp.size() - 1; i >= 0; i--)
- if(cp[i].equal(hh_cp[p])) {isIn = true; break;}
- if(!isIn) cp.push_back(hh_cp[p]);
+ std::pair<DI_Point::Container::iterator,bool> isIn = cp.insert(hh_cp[p]);
+ if(!isIn.second) delete hh_cp[p];
}
+ delete hh;
return signChange;
}
// cut a hex into subtetras along the levelset surface
-void DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Hexa> &Hexas, std::vector<DI_Tetra> &subTetras,
- std::vector<DI_Quad> &surfQuads, std::vector<DI_Triangle> &surfTriangles,
- std::vector<DI_CuttingPoint> &cp, std::vector<DI_QualError> &QError) const
+bool DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Hexa *> &Hexas, std::vector<DI_Tetra *> &subTetras,
+ std::vector<DI_Quad *> &surfQuads, std::vector<DI_Triangle *> &surfTriangles,
+ std::vector<DI_CuttingPoint *> &cp, std::vector<DI_QualError *> &QError)
{
// check if the hex is cut by the level set
int on = 0, pos = 0, neg = 0, ze[8];
for (int i = 0; i < 8; i++){
- if(pt(i).isOnBorder()) ze[on++] = i;
- else if (pt(i).ls() > 0.) pos++;
+ if(pt(i)->isOnBorder()) ze[on++] = i;
+ else if (pt(i)->ls() > 0.) pos++;
else neg++;
}
if (pos && neg) {
- std::vector<DI_Tetra> tetras;
+ std::vector<DI_Tetra *> tetras; tetras.reserve(6);
// Split the quad into sub tetras
splitIntoTetras (tetras);
// each of the sub tetras is split again into sub tetras with respect to the levelset
int nt0 = surfTriangles.size(), nT1, nt1, nT2, nt2;
for (int t = 0; t < (int)tetras.size(); t++) {
nT1 = subTetras.size(); nt1 = surfTriangles.size();
- tetras[t].selfSplit(e, RPNi, subTetras, surfTriangles, cp, QError);
+ tetras[t]->selfSplit(e, RPNi, subTetras, surfTriangles, cp, QError);
nT2 = subTetras.size(); nt2 = surfTriangles.size();
if((nT2 - nT1) == 1 && (nt2 - nt1) == 1){ // only 1 triangle created on a face of a tetra => check if it was not yet created
- if(isLastTrInV(surfTriangles, nt0)) {surfTriangles.pop_back(); nt2--;}
+ if(isLastTrInV(surfTriangles, nt0)) {delete surfTriangles.back(); surfTriangles.pop_back(); nt2--;}
}
}
}
@@ -2813,18 +2881,21 @@ void DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> &RP
if(on == 4){ // the level set is zero on a face of the hex
// assert the nodes are in the same plane
if(!isPlanar(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) {
- printf("Error : The 4 nodes with zero levelset are not planar!\n"); return;}
- // order the 4 nodes
- if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) swap(ze[2], ze[3]);
- // add the quad twice if the face belongs to 2 hexas => remove it later!
- if(ze[0] == 2) surfQuads.push_back(DI_Quad(pt(ze[1]), pt(ze[2]), pt(ze[3]), pt(ze[0]), RPNi.back()->getTag()));
- // to assert that the quad will be split into triangles in the same manner as the hexa into tetras.
- else surfQuads.push_back(DI_Quad(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]), RPNi.back()->getTag()));
+ printf("Error : The 4 nodes with zero levelset are not planar!\n");}
+ else {
+ // order the 4 nodes
+ if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) swap(ze[2], ze[3]);
+ // add the quad twice if the face belongs to 2 hexas => remove it later!
+ if(ze[0] == 2) surfQuads.push_back(new DI_Quad(pt(ze[1]), pt(ze[2]), pt(ze[3]), pt(ze[0]), RPNi.back()->getTag()));
+ // to assert that the quad will be split into triangles in the same manner as the hexa into tetras.
+ else surfQuads.push_back(new DI_Quad(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]), RPNi.back()->getTag()));
+ }
}
for(int i = 0; i < on; i++)
- cp.push_back(DI_CuttingPoint(pt(ze[i])));
- Hexas.push_back(*this);
+ cp.push_back(new DI_CuttingPoint(pt(ze[i])));
+ Hexas.push_back(this);
}
+ return (pos && neg);
}
#endif
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 791fbea19a..f0006bbbd4 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -3,6 +3,7 @@
#include <list>
#include <vector>
+#include <set>
#include <map>
#include <cmath>
#include "DILevelset.h"
@@ -17,6 +18,7 @@
class DI_Point;
class DI_IntegrationPoint;
class DI_CuttingPoint;
+class DI_PointLessThan;
class DI_Element;
class DI_Line; // : public DI_Element
class DI_Triangle; // : public DI_Element
@@ -28,6 +30,7 @@ class DI_QualError;
// --------------------------------------------------------------------------------------------------
class DI_Point
{
+ protected :
// coordinates of the point
double x_, y_, z_;
// vector containing the levelset values of the point
@@ -36,15 +39,15 @@ class DI_Point
// constructors
DI_Point () : x_(0), y_(0), z_(0) {}
DI_Point (double x, double y, double z) : x_(x), y_(y), z_(z) {}
- DI_Point (double x, double y, double z, gLevelset &ls);
+ DI_Point (double x, double y, double z, gLevelset *ls);
DI_Point (double x, double y, double z, const double ls) : x_(x), y_(y), z_(z) {addLs(ls);}
DI_Point (double x, double y, double z, const DI_Element *e,
- const std::vector<const gLevelset *> RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
+ const std::vector<const gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
DI_Point(const DI_Point &p) : x_(p.x()), y_(p.y()), z_(p.z()) {Ls.clear(); Ls = p.Ls;}
// assignment
DI_Point & operator=(const DI_Point & rhs);
// destructor
- ~DI_Point () {Ls.clear();}
+ virtual ~DI_Point () {Ls.clear();}
// add a levelset value (adjusted to 0 if ls<ZERO_LS_TOL) into the vector Ls
void addLs (const double ls);
// add a levelset value evaluated into the element e
@@ -53,15 +56,15 @@ class DI_Point
// delete the last two values and add the chosen one
void chooseLs (const gLevelset *Lsi);
// clear Ls and add the levelset values computed with RPNi
- void computeLs (const DI_Element *e, const std::vector<const gLevelset *> RPNi);
+ void computeLs (const DI_Element *e, const std::vector<const gLevelset *> &RPNi);
// clear Ls and add the levelset value computed with ls
- void computeLs (const gLevelset &ls);
+ void computeLs (const gLevelset *ls);
// remove the last value in Ls and add ls
inline void changeLs (const double ls) {Ls.pop_back(); addLs(ls);}
// change the coordinates
inline void move (double x, double y, double z) {x_ = x; y_ = y; z_ = z;}
// return true if the coordinates of this and p are equal (with a tolerance)
- bool equal(const DI_Point &p) const;
+ bool equal(const DI_Point *p) const;
// return the coordinates
inline double x () const {return x_;}
inline double y () const {return y_;}
@@ -84,12 +87,13 @@ class DI_Point
return 0;
}
// print the coordiantes
- inline void print() const {printf("Point (%g,%g,%g)\n", x_, y_, z_);}
- void printls() const {
+ virtual void print() const {printf("Point (%g,%g,%g)\n", x_, y_, z_);}
+ virtual void printls() const {
printf("Point (%g,%g,%g) ls=(", x_, y_, z_);
for(int i = 0; i < (int)Ls.size(); i++) printf("%g,", Ls[i]);
printf(")\n");
}
+ typedef std::set<DI_Point*,DI_PointLessThan> Container;
};
// --------------------------------------------------------------------------------------------------
class DI_IntegrationPoint
@@ -111,7 +115,7 @@ class DI_IntegrationPoint
// change the value of ls_
inline void setLs (double lsT) {ls_ = lsT;}
// clear Ls and add the levelset values computed with RPNi
- void computeLs (const DI_Element *e, const std::vector<const gLevelset *> RPNi);
+ void computeLs (const DI_Element *e, const std::vector<const gLevelset *> &RPNi);
// change the value of weight_
inline void setWeight (double w) {weight_ = w;}
// return the coordinates
@@ -135,15 +139,43 @@ class DI_IntegrationPoint
printf("IP : x=(%g,%g,%g) xl=(%g,%g,%g) w=%g ls=%g\n", x_, y_, z_, xl_, yl_, zl_, weight_, ls_); }
};
+// --------------------------------------------------------------------------------------------------
+class DI_CuttingPoint : public DI_Point
+{
+ double xl_, yl_, zl_;
+ public:
+ DI_CuttingPoint (const DI_Point *pt);
+ inline void addLocC (double xl, double yl, double zl) {xl_ = xl; yl_ = yl; zl_ = zl;}
+ inline double xl() const {return xl_;}
+ inline double yl() const {return yl_;}
+ inline double zl() const {return zl_;}
+ void print() const {
+ printf("CP : x=(%g,%g,%g) xl=(%g,%g,%g)\n", x_, y_, z_, xl_, yl_, zl_);
+ }
+ void printls() const {
+ printf("CP : x=(%g,%g,%g) xl=(%g,%g,%g) ls=(", x_, y_, z_, xl_, yl_, zl_);
+ for(int i = 0; i < (int)Ls.size(); i++) printf("%g,", Ls[i]);
+ printf(")\n");
+ }
+};
+
+// --------------------------------------------------------------------------------------------------
+class DI_PointLessThan
+{
+ public:
+ static double tolerance;
+ bool operator()(const DI_Point *p1, const DI_Point *p2) const;
+};
+
// --------------------------------------------------------------------------------------------------
// compute the length of a line
static inline double LineLength(double x1, double y1, double z1, double x2, double y2, double z2)
{
return sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2) + (z1 - z2) * (z1 - z2));
}
-static inline double LineLength (const DI_Point &p1, const DI_Point &p2)
+static inline double LineLength (const DI_Point *p1, const DI_Point *p2)
{
- return LineLength(p1.x(), p1.y(), p1.z(), p2.x(), p2.y(), p2.z());
+ return LineLength(p1->x(), p1->y(), p1->z(), p2->x(), p2->y(), p2->z());
}
// compute the surface of a triangle
static inline double TriSurf(double x1, double y1, double z1, double x2, double y2, double z2,
@@ -156,9 +188,9 @@ static inline double TriSurf(double x1, double y1, double z1, double x2, double
+ (z1 * (x2 - x3) - z2 * (x1 - x3) + z3 * (x1 - x2)) *
(z1 * (x2 - x3) - z2 * (x1 - x3) + z3 * (x1 - x2)));
}
-static inline double TriSurf (const DI_Point &p1, const DI_Point &p2, const DI_Point &p3)
+static inline double TriSurf (const DI_Point *p1, const DI_Point *p2, const DI_Point *p3)
{
- return TriSurf(p1.x(), p1.y(), p1.z(), p2.x(), p2.y(), p2.z(), p3.x(), p3.y(), p3.z());
+ return TriSurf(p1->x(), p1->y(), p1->z(), p2->x(), p2->y(), p2->z(), p3->x(), p3->y(), p3->z());
}
// compute the volume of a tetrahedron (base in ccw order to have positive volume)
@@ -171,11 +203,11 @@ static inline double TetraVol(double x1, double y1, double z1, double x2, double
if(vol < 0) {printf("TET HAS NEGATIVE VOLUME = %g\n", vol);}
return vol;
}
-static inline double TetraVol(const DI_Point &p1, const DI_Point &p2,
- const DI_Point &p3, const DI_Point &p4)
+static inline double TetraVol(const DI_Point *p1, const DI_Point *p2,
+ const DI_Point *p3, const DI_Point *p4)
{
- return TetraVol (p1.x(), p1.y(), p1.z(), p2.x(), p2.y(), p2.z(),
- p3.x(), p3.y(), p3.z(), p4.x(), p4.y(), p4.z());
+ return TetraVol (p1->x(), p1->y(), p1->z(), p2->x(), p2->y(), p2->z(),
+ p3->x(), p3->y(), p3->z(), p4->x(), p4->y(), p4->z());
}
// --------------------------------------------------------------------------------------------------
@@ -216,7 +248,7 @@ class DI_Element
int getPolynomialOrder() const {return polOrder_;}
// set the polynomial order of the shape functions
void setPolynomialOrder(int o);
- void setPolynomialOrder(int o, const DI_Element *e, const std::vector<const gLevelset *> RPNi);
+ void setPolynomialOrder(int o, const DI_Element *e, const std::vector<const gLevelset *> &RPNi);
// return tag
int lsTag() const {return lsTag_;}
// return the position of the point with respect to the domain depending on ls_
@@ -238,8 +270,8 @@ class DI_Element
// compute the level set at each vertex of the real element e (same type as this)
// or take i from the map if exist
// and add it to the vertices
- void addLs (const DI_Element *e, const gLevelset &Ls);
- void addLs (const DI_Element *e, const gLevelset &Ls, std::map<int, double> nodeLs[]);
+ void addLs (const DI_Element *e, const gLevelset *Ls);
+ void addLs (const DI_Element *e, const gLevelset *Ls, int iLs, double **nodeLs);
// clear the levelset of the vertices
void clearLs();
// compute the levelset values at the mid edge points and add a quadratic edge,
@@ -247,39 +279,39 @@ class DI_Element
// if xm is too close from the middle of the edgeth edge, do not add the quadratic edge
// if the new quadratic edge create a negative detJ, do not add the quadratic edge and return false
bool addQuadEdge (int edge, DI_Point *xm, const DI_Element *e,
- const std::vector<const gLevelset *> RPNi);
+ const std::vector<const gLevelset *> &RPNi);
// return true if the point pt is inside the element
- bool contain (const DI_Point &pt) const;
+ bool contain (const DI_Point *pt) const;
// return true if the element e is inside the element
// (works only for triangles and quadrangles for the moment)
bool contain (const DI_Element *e) const;
// choose the levelset for each point
void chooseLs (const gLevelset *Lsi);
// map the point pt with local coordinates from the reference element into the element
- DI_Point mappingP (DI_Point &pt) const;
+ void mappingP (DI_Point *pt) const;
// map the DI_IntegrationPoint in with local coordinates from the reference element into the element
- DI_IntegrationPoint mappingIP (DI_IntegrationPoint &in) const;
+ void mappingIP (DI_IntegrationPoint *in) const;
// map the DI_CuttingPoint cp with local coordinates from the reference element into the element
- DI_CuttingPoint mappingCP (DI_CuttingPoint &cp) const;
+ void mappingCP (DI_CuttingPoint *cp) const;
// map the DI_Element el with local coordinates from the reference element into the element
void mappingEl (DI_Element *el) const;
// push into ip the reference integration points to integrate exactly a polynom of order polOrder
virtual void getRefIntegrationPoints (const int polOrder,
- std::vector<DI_IntegrationPoint> &ip) const = 0;
+ std::vector<DI_IntegrationPoint *> &ip) const = 0;
// push into ip the integration points to integrate exactly a polynom of order polOrder
// over the element
- // the local coordiantes of the integration points are computed with loc
+ // the local coordinates of the integration points are computed with loc
// the levelset value at the integration points is computed with e and RPN
void integrationPoints (const int polyOrder, const DI_Element *loc, const DI_Element *e,
- const std::vector<const gLevelset *> RPN,
- std::vector<DI_IntegrationPoint> &ip) const;
+ const std::vector<const gLevelset *> &RPN,
+ std::vector<DI_IntegrationPoint *> &ip) const;
// compute the cutting points on the edges with the last ls in the element
- void getCuttingPoints (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
- std::vector<DI_CuttingPoint> &cp) const;
+ void getCuttingPoints (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_CuttingPoint*> &cp) const;
// return the ith point
- inline DI_Point & pt (int i) const {return i < nbVert() ? *pts_[i] : *mid_[i - nbVert()];}
+ inline DI_Point* pt (int i) const {return i < nbVert() ? pts_[i] : mid_[i - nbVert()];}
// return the ith middle point
- inline DI_Point & mid (int i) const {return *mid_[i];}
+ inline DI_Point* mid (int i) const {return mid_[i];}
// return the coordinates of the ith point
inline double x (int i) const {return i < nbVert() ? pts_[i]->x() : mid_[i - nbVert()]->x();}
inline double y (int i) const {return i < nbVert() ? pts_[i]->y() : mid_[i - nbVert()]->y();}
@@ -308,52 +340,22 @@ class DI_Element
// compute DetJ at each point, return false if 2 values have opposite sign
bool testDetJ() const;
// set the lsTag to +1 if the element is inside the domain (compute in the reference element)
- void computeLsTagDom(const DI_Element *e, const std::vector<const gLevelset *> RPN);
+ void computeLsTagDom(const DI_Element *e, const std::vector<const gLevelset *> &RPN);
// set the lsTag to -1 if the element is not on the boundary of the final levelset
// (compute in the reference element)
- void computeLsTagBound(std::vector<DI_Hexa> &hexas, std::vector<DI_Tetra> &tetras);
- void computeLsTagBound(std::vector<DI_Quad> &quads, std::vector<DI_Triangle> &triangles);
+ void computeLsTagBound(std::vector<DI_Hexa *> &hexas, std::vector<DI_Tetra *> &tetras);
+ void computeLsTagBound(std::vector<DI_Quad *> &quads, std::vector<DI_Triangle *> &triangles);
// print the coordinates of the points of the element
void print () const;
// print the coordinates and the levelset values of the points of the element
void printls () const;
};
-// --------------------------------------------------------------------------------------------------
-class DI_CuttingPoint
+class DI_ElementLessThan
{
- double x_, y_, z_;
- double xl_, yl_, zl_;
- std::vector<double> Ls;
public:
- DI_CuttingPoint (const DI_Point &pt);
- inline void addLocC (double xl, double yl, double zl) {xl_ = xl; yl_ = yl; zl_ = zl;}
- inline void move (double x, double y, double z) {x_ = x; y_ = y; z_ = z;}
- void addLs (const double ls);
- void addLs (const DI_Element *e);
- inline double x () const {return x_;}
- inline double y () const {return y_;}
- inline double z () const {return z_;}
- inline double xl() const {return xl_;}
- inline double yl() const {return yl_;}
- inline double zl() const {return zl_;}
- inline double ls() const {return Ls.back();}
- inline double ls(int i) const {return Ls[i];}
- // return the position of the point with respect to the domain depending on the last value in Ls
- inline bool isInsideDomain () const {return Ls.back() < 0.;}
- inline bool isOutsideDomain () const {return Ls.back() > 0.;}
- inline bool isOnBorder () const {return Ls.back() == 0.;}
- inline int sizeLs() const {return Ls.size();}
- void chooseLs (const gLevelset *Lsi);
- bool equal (const DI_CuttingPoint &p) const;
- inline void print() const {
- printf("CP : x=(%g,%g,%g) xl=(%g,%g,%g)\n", x_, y_, z_, xl_, yl_, zl_);
- }
- void printls() const {
- printf("CP : x=(%g,%g,%g) xl=(%g,%g,%g) ls=(", x_, y_, z_, xl_, yl_, zl_);
- for(int i = 0; i < (int)Ls.size(); i++) printf("%g,", Ls[i]);
- printf(")\n");
- }
+ static double tolerance;
+ bool operator()(const DI_Element *v1, const DI_Element *v2) const;
};
// --------------------------------------------------------------------------------------------------
@@ -380,12 +382,12 @@ class DI_Line : public DI_Element
pts_[1] = new DI_Point(x1, y1, z1);
integral_ = length;
}
- DI_Line (const DI_Point &pt0, const DI_Point &pt1, const int tag = -1)
+ DI_Line (const DI_Point *pt0, const DI_Point *pt1, const int tag = -1)
{
lsTag_ = tag;
pts_ = new DI_Point*[2];
- pts_[0] = new DI_Point(pt0);
- pts_[1] = new DI_Point(pt1);
+ pts_[0] = new DI_Point(*pt0);
+ pts_[1] = new DI_Point(*pt1);
integral_ = LineLength(pt0, pt1);
}
~DI_Line () {
@@ -407,7 +409,7 @@ class DI_Line : public DI_Element
virtual void computeIntegral();
virtual double refIntegral() const {return 2.;}
virtual void getRefIntegrationPoints (const int polynomialOrder,
- std::vector<DI_IntegrationPoint> &ipS) const;
+ std::vector<DI_IntegrationPoint *> &ipS) const;
inline void vert(const int edge, int &s1, int &s2) const {
s1 = 0; s2 = 1;}
void midV (const int e, int *s, int &n) const {
@@ -417,13 +419,13 @@ class DI_Line : public DI_Element
void getGradShapeFunctions (const double u, const double v, const double w,
double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
- bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_CuttingPoint> &cp, const int polynomialOrderL,
- std::vector<DI_Line> &subLines, int recurLevel = 0, std::map<int, double> nodeLs[2] = NULL) const;
- void cut(const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cp) const;
+ bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ DI_Point::Container &cp, const int polynomialOrderL,
+ std::vector<DI_Line *> &subLines, int recurLevel = 0, double **nodeLs = NULL) const;
+ bool cut(const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint*> &cp);
void selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cuttingPoints) const;
+ std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint*> &cuttingPoints) const;
inline double length() const {return integral_;}
};
@@ -457,13 +459,13 @@ class DI_Triangle : public DI_Element
pts_[2] = new DI_Point(x2, y2, z2);
integral_ = surface;
}
- DI_Triangle (const DI_Point &pt0, const DI_Point &pt1, const DI_Point &pt2, const int tag = -1)
+ DI_Triangle (const DI_Point *pt0, const DI_Point *pt1, const DI_Point *pt2, const int tag = -1)
{
lsTag_ = tag;
pts_ = new DI_Point*[3];
- pts_[0] = new DI_Point(pt0);
- pts_[1] = new DI_Point(pt1);
- pts_[2] = new DI_Point(pt2);
+ pts_[0] = new DI_Point(*pt0);
+ pts_[1] = new DI_Point(*pt1);
+ pts_[2] = new DI_Point(*pt2);
integral_ = TriSurf(pt0,pt1,pt2);
}
~DI_Triangle () {
@@ -485,7 +487,7 @@ class DI_Triangle : public DI_Element
virtual void computeIntegral();
virtual double refIntegral() const {return 0.5;}
virtual void getRefIntegrationPoints (const int polynomialOrder,
- std::vector<DI_IntegrationPoint> &ipS) const;
+ std::vector<DI_IntegrationPoint *> &ipS) const;
inline void vert(const int edge, int &s1, int &s2) const {
int v[3][2] = {{0, 1}, {1, 2}, {2, 0}};
s1 = v[edge][0]; s2 = v[edge][1];
@@ -502,18 +504,18 @@ class DI_Triangle : public DI_Element
void getGradShapeFunctions (const double u, const double v, const double w,
double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
- bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+ bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, int recurLevel = 0, std::map<int, double> nodeLs[3] = NULL) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cp) const;
- void splitIntoSubTriangles (std::vector<DI_Triangle> &triangles) const;
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, int recurLevel = 0, double **nodeLs = NULL) const;
+ bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint*> &cp);
+ void splitIntoSubTriangles (std::vector<DI_Triangle *> &triangles) const;
void selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cuttingPoints) const;
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint*> &cuttingPoints) const;
double quality () const;
inline double surface() const {return integral_;}
};
@@ -551,15 +553,15 @@ class DI_Quad : public DI_Element
pts_[3] = new DI_Point(x3, y3, z3);
integral_ = surf;
}
- DI_Quad (const DI_Point &pt0, const DI_Point &pt1, const DI_Point &pt2, const DI_Point &pt3,
+ DI_Quad (const DI_Point *pt0, const DI_Point *pt1, const DI_Point *pt2, const DI_Point *pt3,
const int tag = -1)
{
lsTag_ = tag;
pts_ = new DI_Point*[4];
- pts_[0] = new DI_Point(pt0);
- pts_[1] = new DI_Point(pt1);
- pts_[2] = new DI_Point(pt2);
- pts_[3] = new DI_Point(pt3);
+ pts_[0] = new DI_Point(*pt0);
+ pts_[1] = new DI_Point(*pt1);
+ pts_[2] = new DI_Point(*pt2);
+ pts_[3] = new DI_Point(*pt3);
integral_ = TriSurf(pt0, pt1, pt2) + TriSurf(pt0, pt2, pt3);
}
~DI_Quad () {
@@ -581,7 +583,7 @@ class DI_Quad : public DI_Element
virtual void computeIntegral();
virtual double refIntegral() const {return 4.;}
virtual void getRefIntegrationPoints (const int polynomialOrder,
- std::vector<DI_IntegrationPoint> &ipS) const;
+ std::vector<DI_IntegrationPoint *> &ipS) const;
inline void vert(const int edge, int &s1, int &s2) const{
int v[4][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}};
s1 = v[edge][0]; s2 = v[edge][1];
@@ -600,15 +602,15 @@ class DI_Quad : public DI_Element
void getGradShapeFunctions (const double u, const double v, const double w,
double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
- bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+ bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, int recurLevel = 0, std::map<int, double> nodeLs[4] = NULL) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
- std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cp) const;
- void splitIntoTriangles (std::vector<DI_Triangle> &triangles) const;
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, int recurLevel = 0, double **nodeLs = NULL) const;
+ bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
+ std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint*> &cp);
+ void splitIntoTriangles (std::vector<DI_Triangle *> &triangles) const;
inline double surface() const {return integral_;}
};
@@ -647,13 +649,13 @@ class DI_Tetra : public DI_Element
pts_[3] = new DI_Point(x3, y3, z3);
integral_ = vol;
}
- DI_Tetra (const DI_Point &pt0, const DI_Point &pt1, const DI_Point &pt2, const DI_Point &pt3)
+ DI_Tetra (const DI_Point *pt0, const DI_Point *pt1, const DI_Point *pt2, const DI_Point *pt3)
{
pts_ = new DI_Point*[4];
- pts_[0] = new DI_Point(pt0);
- pts_[1] = new DI_Point(pt1);
- pts_[2] = new DI_Point(pt2);
- pts_[3] = new DI_Point(pt3);
+ pts_[0] = new DI_Point(*pt0);
+ pts_[1] = new DI_Point(*pt1);
+ pts_[2] = new DI_Point(*pt2);
+ pts_[3] = new DI_Point(*pt3);
integral_ = TetraVol(pt0, pt1, pt2, pt3);
}
~DI_Tetra () {
@@ -675,7 +677,7 @@ class DI_Tetra : public DI_Element
virtual void computeIntegral();
virtual double refIntegral() const {return 1. / 6.;}
virtual void getRefIntegrationPoints (const int polynomialOrder,
- std::vector<DI_IntegrationPoint> &ip) const;
+ std::vector<DI_IntegrationPoint *> &ip) const;
inline void vert(const int edge, int &s1, int &s2) const {
int v[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
s1 = v[edge][0]; s2 = v[edge][1];
@@ -695,18 +697,18 @@ class DI_Tetra : public DI_Element
void getGradShapeFunctions (const double u, const double v, const double w,
double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
- bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+ bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderT, const int polynomialOrderQ, const int polynomialOrderTr,
- std::vector<DI_Tetra> &subTetras, std::vector<DI_Quad> &surfQuads,
- std::vector<DI_Triangle> &surfTriangles, int recurLevel = 0, std::map<int, double> nodeLs[4] = NULL) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Tetra> &subTetras, std::vector<DI_Quad> &surfQuads,
- std::vector<DI_Triangle> &surfTriangles, std::vector<DI_CuttingPoint> &cp,
- std::vector<DI_QualError> &QE) const;
+ std::vector<DI_Tetra *> &subTetras, std::vector<DI_Quad *> &surfQuads,
+ std::vector<DI_Triangle *> &surfTriangles, int recurLevel = 0, double **nodeLs = NULL) const;
+ bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Tetra *> &subTetras, std::vector<DI_Quad *> &surfQuads,
+ std::vector<DI_Triangle *> &surfTriangles, std::vector<DI_CuttingPoint*> &cp,
+ std::vector<DI_QualError *> &QE);
void selfSplit ( const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Tetra> &subTetras, std::vector<DI_Triangle> &surfTriangles,
- std::vector<DI_CuttingPoint> &cuttingPoints, std::vector<DI_QualError> &QE) const;
+ std::vector<DI_Tetra *> &subTetras, std::vector<DI_Triangle *> &surfTriangles,
+ std::vector<DI_CuttingPoint*> &cuttingPoints, std::vector<DI_QualError *> &QE) const;
double quality () const;
inline double volume() const {return integral_;}
};
@@ -759,13 +761,13 @@ class DI_Hexa : public DI_Element
pts_[6] = new DI_Point(x6, y6, z6); pts_[7] = new DI_Point(x7, y7, z7);
integral_ = vol;
}
- DI_Hexa (const DI_Point &pt0, const DI_Point &pt1, const DI_Point &pt2, const DI_Point &pt3,
- const DI_Point &pt4, const DI_Point &pt5, const DI_Point &pt6, const DI_Point &pt7) {
+ DI_Hexa (const DI_Point *pt0, const DI_Point *pt1, const DI_Point *pt2, const DI_Point *pt3,
+ const DI_Point *pt4, const DI_Point *pt5, const DI_Point *pt6, const DI_Point *pt7) {
pts_ = new DI_Point*[8];
- pts_[0] = new DI_Point(pt0); pts_[1] = new DI_Point(pt1);
- pts_[2] = new DI_Point(pt2); pts_[3] = new DI_Point(pt3);
- pts_[4] = new DI_Point(pt4); pts_[5] = new DI_Point(pt5);
- pts_[6] = new DI_Point(pt6); pts_[7] = new DI_Point(pt7);
+ pts_[0] = new DI_Point(*pt0); pts_[1] = new DI_Point(*pt1);
+ pts_[2] = new DI_Point(*pt2); pts_[3] = new DI_Point(*pt3);
+ pts_[4] = new DI_Point(*pt4); pts_[5] = new DI_Point(*pt5);
+ pts_[6] = new DI_Point(*pt6); pts_[7] = new DI_Point(*pt7);
integral_ = TetraVol(pt0, pt1, pt3, pt4) + TetraVol(pt1, pt4, pt5, pt7) +
TetraVol(pt1, pt3, pt4, pt7) + TetraVol(pt2, pt5, pt6, pt7) +
TetraVol(pt1, pt2, pt3, pt7) + TetraVol(pt1, pt5, pt2, pt7);
@@ -789,7 +791,7 @@ class DI_Hexa : public DI_Element
virtual void computeIntegral();
virtual double refIntegral() const {return 8.;}
virtual void getRefIntegrationPoints (const int polynomialOrder,
- std::vector<DI_IntegrationPoint> &ip) const;
+ std::vector<DI_IntegrationPoint *> &ip) const;
inline void vert(const int edge, int &s1, int &s2) const {
int v[12][2] = {{0, 1}, {1, 2}, {2, 3}, {3, 0}, {0, 4}, {1, 5},
{2, 6}, {3, 7}, {4, 5}, {5, 6}, {6, 7}, {7, 4}};
@@ -824,18 +826,18 @@ class DI_Hexa : public DI_Element
void getGradShapeFunctions (const double u, const double v, const double w,
double s[][3], int order = -1) const;
double detJ (const double &xP, const double &yP, const double &zP) const;
- bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
- std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
+ bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ std::vector<DI_IntegrationPoint *> &ipS, DI_Point::Container &cp,
const int polynomialOrderH, const int polynomialOrderT,
const int polynomialOrderQ, const int polynomialOrderTr,
- std::vector<DI_Hexa> ¬CutHexas, std::vector<DI_Tetra> &subTetras,
- std::vector<DI_Quad> &surfQuads, std::vector<DI_Triangle> &surfTriangles,
- std::vector<DI_Line> &frontLines, int recurLevel = 0, std::map<int, double> nodeLs[8] = NULL) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
- std::vector<DI_Hexa> &unCutHexas, std::vector<DI_Tetra> &subTetras,
- std::vector<DI_Quad> &surfQuads, std::vector<DI_Triangle> &surfTriangles,
- std::vector<DI_CuttingPoint> &cp, std::vector<DI_QualError> &QE) const;
- void splitIntoTetras(std::vector<DI_Tetra> &tetras) const;
+ std::vector<DI_Hexa *> ¬CutHexas, std::vector<DI_Tetra *> &subTetras,
+ std::vector<DI_Quad *> &surfQuads, std::vector<DI_Triangle *> &surfTriangles,
+ std::vector<DI_Line *> &frontLines, int recurLevel = 0, double **nodeLs = NULL) const;
+ bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ std::vector<DI_Hexa *> &unCutHexas, std::vector<DI_Tetra *> &subTetras,
+ std::vector<DI_Quad *> &surfQuads, std::vector<DI_Triangle *> &surfTriangles,
+ std::vector<DI_CuttingPoint*> &cp, std::vector<DI_QualError *> &QE);
+ void splitIntoTetras(std::vector<DI_Tetra *> &tetras) const;
inline double volume() const {return integral_;}
};
@@ -845,21 +847,31 @@ class DI_Hexa : public DI_Element
// -------------------------------------------------------------------------------------------------
class DI_QualError
{
- DI_Point p11_, p12_, p21_, p22_;
+ DI_Point **pts_;
public:
- DI_QualError (DI_Point p11, DI_Point p12, DI_Point p21, DI_Point p22)
- : p11_(p11), p12_(p12), p21_(p21), p22_(p22) {}
- inline DI_Point pt (int i) const {
- if(i == 0) return p11_;
- if(i == 1) return p21_;
- if(i == 2) return p12_;
- if(i == 3) return p22_;
+ DI_QualError (DI_Point *p11, DI_Point *p12, DI_Point *p21, DI_Point *p22) {
+ pts_ = new DI_Point*[4];
+ pts_[0] = new DI_Point(*p11);
+ pts_[1] = new DI_Point(*p12);
+ pts_[2] = new DI_Point(*p21);
+ pts_[3] = new DI_Point(*p22);
+ }
+ ~DI_QualError () {
+ for(int i = 0; i < 4; i++)
+ delete pts_[i];
+ if(pts_) delete[] pts_;
+ }
+ inline DI_Point* pt (int i) const {
+ if(i == 0) return pts_[0];
+ if(i == 1) return pts_[1];
+ if(i == 2) return pts_[2];
+ if(i == 3) return pts_[3];
printf("DI_QualError::pt only accept indices from 0 to 3!\n");
- DI_Point p; return p;
+ DI_Point* p = NULL; return p;
}
void print(const DI_Element *e) const{
- DI_Point pt1 = p11_, pt2 = p12_, pt3 = p21_, pt4 = p22_;
- e->mappingP(pt1); e->mappingP(pt2); e->mappingP(pt3); e->mappingP(pt4);
+ DI_Point pt1(*pts_[0]), pt2(*pts_[1]), pt3(*pts_[2]), pt4(*pts_[3]);
+ e->mappingP(&pt1); e->mappingP(&pt2); e->mappingP(&pt3); e->mappingP(&pt4);
printf("Cannot assert best quality for the last face of the Prism \n");
printf("=> edges (%g,%g,%g),(%g,%g,%g) and (%g,%g,%g),(%g,%g,%g) may cross in ",
pt1.x(), pt1.y(), pt1.z(), pt2.x(), pt2.y(), pt2.z(),
diff --git a/contrib/DiscreteIntegration/recurCut.cpp b/contrib/DiscreteIntegration/recurCut.cpp
index a54f06376d..aeb6a80c4e 100644
--- a/contrib/DiscreteIntegration/recurCut.cpp
+++ b/contrib/DiscreteIntegration/recurCut.cpp
@@ -3,8 +3,8 @@
#include "recurCut.h"
-inline DI_Point middle (const DI_Point &p1, const DI_Point &p2) {
- return DI_Point ((p1.x() + p2.x()) / 2, (p1.y() + p2.y()) / 2, (p1.z() + p2.z()) / 2);
+inline DI_Point* middle (const DI_Point *p1, const DI_Point *p2) {
+ return new DI_Point ((p1->x() + p2->x()) / 2, (p1->y() + p2->y()) / 2, (p1->z() + p2->z()) / 2);
}
// create sub-elements recursively (according to a pattern) until maxlevel
@@ -14,187 +14,241 @@ void recurCut(RecurElement *re, int maxlevel, int level)
if(re->type() == DI_LIN) {
// p1 p12 p2
- DI_Point p1 = re->el->pt(0);
- DI_Point p2 = re->el->pt(1);
- DI_Point p12 = middle(p1, p2);
- RecurElement *re0 = new RecurElement(&DI_Line(p1, p12));
+ DI_Point* p1 = re->el->pt(0);
+ DI_Point* p2 = re->el->pt(1);
+ DI_Point* p12 = middle(p1, p2);
+ DI_Line *l0 = new DI_Line(p1, p12);
+ RecurElement *re0 = new RecurElement(l0);
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Line(p12, p2));
+ DI_Line *l1 = new DI_Line(p12, p2);
+ RecurElement *re1 = new RecurElement(l1);
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
+ delete p12;
+ delete l0; delete l1;
}
if(re->type() == DI_TRI) {
// p3
// p13 p23
// p1 p12 p2
- DI_Point p1 = re->el->pt(0);
- DI_Point p2 = re->el->pt(1);
- DI_Point p3 = re->el->pt(2);
- DI_Point p12 = middle(p1, p2);
- DI_Point p13 = middle(p1, p3);
- DI_Point p23 = middle(p2, p3);
- RecurElement *re0 = new RecurElement(&DI_Triangle(p1, p12, p13));
+ DI_Point *p1 = re->el->pt(0);
+ DI_Point *p2 = re->el->pt(1);
+ DI_Point *p3 = re->el->pt(2);
+ DI_Point *p12 = middle(p1, p2);
+ DI_Point *p13 = middle(p1, p3);
+ DI_Point *p23 = middle(p2, p3);
+ DI_Triangle *t0 = new DI_Triangle(p1, p12, p13);
+ RecurElement *re0 = new RecurElement(t0);
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Triangle(p2, p23, p12));
+ DI_Triangle *t1 = new DI_Triangle(p2, p23, p12);
+ RecurElement *re1 = new RecurElement(t1);
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
- RecurElement *re2 = new RecurElement(&DI_Triangle(p3, p13, p23));
+ DI_Triangle *t2 = new DI_Triangle(p3, p13, p23);
+ RecurElement *re2 = new RecurElement(t2);
recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re;
- RecurElement *re3 = new RecurElement(&DI_Triangle(p12, p23, p13));
+ DI_Triangle *t3 = new DI_Triangle(p12, p23, p13);
+ RecurElement *re3 = new RecurElement(t3);
recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re;
+ delete p12; delete p13; delete p23;
+ delete t0; delete t1; delete t2; delete t3;
}
else if(re->type() == DI_QUA) {
// p4 p34 p3
// p14 p1234 p23
// p1 p12 p2
- DI_Point p1 = re->el->pt(0);
- DI_Point p2 = re->el->pt(1);
- DI_Point p3 = re->el->pt(2);
- DI_Point p4 = re->el->pt(3);
- DI_Point p12 = middle(p1, p2);
- DI_Point p23 = middle(p2, p3);
- DI_Point p34 = middle(p3, p4);
- DI_Point p14 = middle(p1, p4);
- DI_Point p1234 = middle(p12, p34);
- RecurElement *re0 = new RecurElement(&DI_Quad(p1, p12, p1234, p14));
+ DI_Point *p1 = re->el->pt(0);
+ DI_Point *p2 = re->el->pt(1);
+ DI_Point *p3 = re->el->pt(2);
+ DI_Point *p4 = re->el->pt(3);
+ DI_Point *p12 = middle(p1, p2);
+ DI_Point *p23 = middle(p2, p3);
+ DI_Point *p34 = middle(p3, p4);
+ DI_Point *p14 = middle(p1, p4);
+ DI_Point *p1234 = middle(p12, p34);
+ DI_Quad *q0 = new DI_Quad(p1, p12, p1234, p14);
+ RecurElement *re0 = new RecurElement(q0);
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Quad(p12, p2, p23, p1234));
+ DI_Quad *q1 = new DI_Quad(p12, p2, p23, p1234);
+ RecurElement *re1 = new RecurElement(q1);
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
- RecurElement *re2 = new RecurElement(&DI_Quad(p1234, p23, p3, p34));
+ DI_Quad *q2 = new DI_Quad(p1234, p23, p3, p34);
+ RecurElement *re2 = new RecurElement(q2);
recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re;
- RecurElement *re3 = new RecurElement(&DI_Quad(p14, p1234, p34, p4));
+ DI_Quad *q3 = new DI_Quad(p14, p1234, p34, p4);
+ RecurElement *re3 = new RecurElement(q3);
recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re;
+ delete p12; delete p23; delete p34;
+ delete p14; delete p1234;
+ delete q0; delete q1; delete q2; delete q3;
}
else if(re->type() == DI_TET) {
- DI_Point p1 = re->el->pt(0);
- DI_Point p2 = re->el->pt(1);
- DI_Point p3 = re->el->pt(2);
- DI_Point p4 = re->el->pt(3);
- DI_Point p12 = middle(p1, p2);
- DI_Point p13 = middle(p1, p3);
- DI_Point p14 = middle(p1, p4);
- DI_Point p23 = middle(p2, p3);
- DI_Point p24 = middle(p2, p4);
- DI_Point p34 = middle(p3, p4);
- RecurElement *re0 = new RecurElement(&DI_Tetra(p1, p12, p13, p14));
+ // p4
+ // p14 p34
+ // p24 p1 p13 p3
+ // p12 p23
+ // p2
+ DI_Point *p1 = re->el->pt(0);
+ DI_Point *p2 = re->el->pt(1);
+ DI_Point *p3 = re->el->pt(2);
+ DI_Point *p4 = re->el->pt(3);
+ DI_Point *p12 = middle(p1, p2);
+ DI_Point *p13 = middle(p1, p3);
+ DI_Point *p14 = middle(p1, p4);
+ DI_Point *p23 = middle(p2, p3);
+ DI_Point *p24 = middle(p2, p4);
+ DI_Point *p34 = middle(p3, p4);
+ DI_Tetra *t0 = new DI_Tetra(p1, p12, p13, p14);
+ RecurElement *re0 = new RecurElement(t0);
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Tetra(p2, p23, p12, p24));
+ DI_Tetra *t1 = new DI_Tetra(p2, p23, p12, p24);
+ RecurElement *re1 = new RecurElement(t1);
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
- RecurElement *re2 = new RecurElement(&DI_Tetra(p3, p13, p23, p34));
+ DI_Tetra *t2 = new DI_Tetra(p3, p13, p23, p34);
+ RecurElement *re2 = new RecurElement(t2);
recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re;
- RecurElement *re3 = new RecurElement(&DI_Tetra(p4, p14, p34, p24));
+ DI_Tetra *t3 = new DI_Tetra(p4, p14, p34, p24);
+ RecurElement *re3 = new RecurElement(t3);
recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re;
- RecurElement *re4 = new RecurElement(&DI_Tetra(p12, p14, p24, p34));
+ DI_Tetra *t4 = new DI_Tetra(p12, p14, p24, p34);
+ RecurElement *re4 = new RecurElement(t4);
recurCut(re4, maxlevel, level);
re->sub[4] = re4; re4->super = re;
- RecurElement *re5 = new RecurElement(&DI_Tetra(p12, p23, p34, p24));
+ DI_Tetra *t5 = new DI_Tetra(p12, p23, p34, p24);
+ RecurElement *re5 = new RecurElement(t5);
recurCut(re5, maxlevel, level);
re->sub[5] = re5; re5->super = re;
- RecurElement *re6 = new RecurElement(&DI_Tetra(p12, p13, p34, p23));
+ DI_Tetra *t6 = new DI_Tetra(p12, p13, p34, p23);
+ RecurElement *re6 = new RecurElement(t6);
recurCut(re6, maxlevel, level);
re->sub[6] = re6; re6->super = re;
- RecurElement *re7 = new RecurElement(&DI_Tetra(p12, p13, p14, p34));
+ DI_Tetra *t7 = new DI_Tetra(p12, p13, p14, p34);
+ RecurElement *re7 = new RecurElement(t7);
recurCut(re7, maxlevel, level);
re->sub[7] = re7; re7->super = re;
+ delete p12; delete p13; delete p14;
+ delete p23; delete p24; delete p34;
+ delete t0; delete t1; delete t2; delete t3;
+ delete t4; delete t5; delete t6; delete t7;
}
else if(re->type() == DI_HEX) {
- DI_Point p1 = re->el->pt(0);
- DI_Point p2 = re->el->pt(1);
- DI_Point p3 = re->el->pt(2);
- DI_Point p4 = re->el->pt(3);
- DI_Point p5 = re->el->pt(4);
- DI_Point p6 = re->el->pt(5);
- DI_Point p7 = re->el->pt(6);
- DI_Point p8 = re->el->pt(7);
- DI_Point p12 = middle(p1, p2);
- DI_Point p14 = middle(p1, p4);
- DI_Point p15 = middle(p1, p5);
- DI_Point p23 = middle(p2, p3);
- DI_Point p26 = middle(p2, p6);
- DI_Point p34 = middle(p3, p4);
- DI_Point p37 = middle(p3, p7);
- DI_Point p48 = middle(p4, p8);
- DI_Point p56 = middle(p5, p6);
- DI_Point p58 = middle(p5, p8);
- DI_Point p67 = middle(p6, p7);
- DI_Point p78 = middle(p7, p8);
- DI_Point p1234 = middle(p12, p34);
- DI_Point p1256 = middle(p12, p56);
- DI_Point p1458 = middle(p14, p58);
- DI_Point p2367 = middle(p23, p67);
- DI_Point p3478 = middle(p34, p78);
- DI_Point p5678 = middle(p56, p78);
- DI_Point p12345678 = middle(p1234, p5678);
- RecurElement *re0 = new RecurElement(&DI_Hexa(p1, p12, p1234, p14, p15, p1256, p12345678, p1458));
+ DI_Point *p1 = re->el->pt(0);
+ DI_Point *p2 = re->el->pt(1);
+ DI_Point *p3 = re->el->pt(2);
+ DI_Point *p4 = re->el->pt(3);
+ DI_Point *p5 = re->el->pt(4);
+ DI_Point *p6 = re->el->pt(5);
+ DI_Point *p7 = re->el->pt(6);
+ DI_Point *p8 = re->el->pt(7);
+ DI_Point *p12 = middle(p1, p2);
+ DI_Point *p14 = middle(p1, p4);
+ DI_Point *p15 = middle(p1, p5);
+ DI_Point *p23 = middle(p2, p3);
+ DI_Point *p26 = middle(p2, p6);
+ DI_Point *p34 = middle(p3, p4);
+ DI_Point *p37 = middle(p3, p7);
+ DI_Point *p48 = middle(p4, p8);
+ DI_Point *p56 = middle(p5, p6);
+ DI_Point *p58 = middle(p5, p8);
+ DI_Point *p67 = middle(p6, p7);
+ DI_Point *p78 = middle(p7, p8);
+ DI_Point *p1234 = middle(p12, p34);
+ DI_Point *p1256 = middle(p12, p56);
+ DI_Point *p1458 = middle(p14, p58);
+ DI_Point *p2367 = middle(p23, p67);
+ DI_Point *p3478 = middle(p34, p78);
+ DI_Point *p5678 = middle(p56, p78);
+ DI_Point *p12345678 = middle(p1234, p5678);
+ DI_Hexa *h0 = new DI_Hexa(p1, p12, p1234, p14, p15, p1256, p12345678, p1458);
+ RecurElement *re0 = new RecurElement(h0);
recurCut(re0, maxlevel, level);
re->sub[0] = re0; re0->super = re;
- RecurElement *re1 = new RecurElement(&DI_Hexa(p12, p2, p23, p1234, p1256, p26, p2367, p12345678));
+ DI_Hexa *h1 = new DI_Hexa(p12, p2, p23, p1234, p1256, p26, p2367, p12345678);
+ RecurElement *re1 = new RecurElement(h1);
recurCut(re1, maxlevel, level);
re->sub[1] = re1; re1->super = re;
- RecurElement *re2 = new RecurElement(&DI_Hexa(p1234, p23, p3, p34, p12345678, p2367, p37, p3478));
+ DI_Hexa *h2 = new DI_Hexa(p1234, p23, p3, p34, p12345678, p2367, p37, p3478);
+ RecurElement *re2 = new RecurElement(h2);
recurCut(re2, maxlevel, level);
re->sub[2] = re2; re2->super = re;
- RecurElement *re3 = new RecurElement(&DI_Hexa(p14, p1234, p34, p4, p1458, p12345678, p3478, p48));
+ DI_Hexa *h3 = new DI_Hexa(p14, p1234, p34, p4, p1458, p12345678, p3478, p48);
+ RecurElement *re3 = new RecurElement(h3);
recurCut(re3, maxlevel, level);
re->sub[3] = re3; re3->super = re;
- RecurElement *re4 = new RecurElement(&DI_Hexa(p15, p1256, p12345678, p1458, p5, p56, p5678, p58));
+ DI_Hexa *h4 = new DI_Hexa(p15, p1256, p12345678, p1458, p5, p56, p5678, p58);
+ RecurElement *re4 = new RecurElement(h4);
recurCut(re4, maxlevel, level);
re->sub[4] = re4; re4->super = re;
- RecurElement *re5 = new RecurElement(&DI_Hexa(p1256, p26, p2367, p12345678, p56, p6, p67, p5678));
+ DI_Hexa *h5 = new DI_Hexa(p1256, p26, p2367, p12345678, p56, p6, p67, p5678);
+ RecurElement *re5 = new RecurElement(h5);
recurCut(re5, maxlevel, level);
re->sub[5] = re5; re5->super = re;
- RecurElement *re6 = new RecurElement(&DI_Hexa(p12345678, p2367, p37, p3478, p5678, p67, p7, p78));
+ DI_Hexa *h6 = new DI_Hexa(p12345678, p2367, p37, p3478, p5678, p67, p7, p78);
+ RecurElement *re6 = new RecurElement(h6);
recurCut(re6, maxlevel, level);
re->sub[6] = re6; re6->super = re;
- RecurElement *re7 = new RecurElement(&DI_Hexa(p1458, p12345678, p3478, p48, p58, p5678, p78, p8));
+ DI_Hexa *h7 = new DI_Hexa(p1458, p12345678, p3478, p48, p58, p5678, p78, p8);
+ RecurElement *re7 = new RecurElement(h7);
recurCut(re7, maxlevel, level);
re->sub[7] = re7; re7->super = re;
+ delete p12; delete p14; delete p15;
+ delete p23; delete p26; delete p34;
+ delete p37; delete p48; delete p56;
+ delete p58; delete p67; delete p78;
+ delete p1234; delete p1256;
+ delete p1458; delete p2367;
+ delete p3478; delete p5678;
+ delete p12345678;
+ delete h0; delete h1; delete h2; delete h3;
+ delete h4; delete h5; delete h6; delete h7;
}
}
// return true if the element re->el is crossed or run along by a primitive levelset in RPN and by the final levelset
// (the levelset is computed with the values at the nodes of the element e)
bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const gLevelset *> &RPN,
- std::map<int, double> nodeLs[8]) {
+ double **nodeLs) {
bool cS = false;
DI_Element* elem = re->root()->el;
- std::vector<DI_CuttingPoint> cp;
+ std::vector<DI_CuttingPoint *> cp;
std::vector<const gLevelset *> RPNi;
for(unsigned int l = 0; l < RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- elem->addLs(e, *Lsi, nodeLs);
+ elem->addLs(e, Lsi, l, nodeLs);
for(unsigned int i = 0; i < cp.size(); i++)
- cp[i].addLs(elem);
+ cp[i]->addLs(elem);
if (re->super) re->el->addLs(elem);
re->el->getCuttingPoints(elem, RPNi, cp);
}
else {
for(unsigned int p = 0; p < cp.size(); p++)
- cp[p].chooseLs(Lsi);
+ cp[p]->chooseLs(Lsi);
if (re->super) re->el->chooseLs(Lsi);
}
}
- re->el->clearLs();
re->root()->el->clearLs();
+ if(re->super) re->el->clearLs();
if(cp.size() > 1 || (re->type() == DI_LIN && cp.size() > 0)) {
for(int i = 0; i < (int)cp.size(); i++)
- if(cp[i].ls() == 0) {cS = true; break;}
+ if(cp[i]->ls() == 0) {cS = true; break;}
}
+ for(int i = 0; i < (int)cp.size(); i++)
+ delete cp[i];
return cS;
}
@@ -202,7 +256,7 @@ bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const
// If it has no sub RecurElement, set isCrossed to true if the element is crossed or run along by the levelset
//(the levelset is computed with the values at the nodes of the triangle e)
bool computeIsCrossed (RecurElement *re, const DI_Element *e, const std::vector<const gLevelset *> &RPN,
- std::map<int, double> nodeLs[8]) {
+ double **nodeLs) {
if (!re->sub[0])
re->isCrossed = signChange(re, e, RPN, nodeLs);
else {
@@ -236,20 +290,18 @@ void recurChangeVisibility(RecurElement *re, const std::vector<const gLevelset *
}
else {
re->el->printls();
- double v = re->Ls();
+ double v = re->ls();
double vm = 0.;
if(!re->sub[0]->sub[0]) {
for(int i = 0; i < re->nbSub(); i++){
- re->sub[i]->el->printls();
- vm += re->sub[i]->Ls();
+ vm += re->sub[i]->ls();
}
vm = vm / re->nbSub();
}
else {
for(int i = 0; i < re->nbSub(); i++){
for(int j = 0; j < re->sub[0]->nbSub(); j++){
- re->sub[i]->sub[j]->el->printls();
- vm += re->sub[i]->sub[j]->Ls();
+ vm += re->sub[i]->sub[j]->ls();
}
}
vm = vm / (re->nbSub()*re->sub[0]->nbSub());
@@ -260,7 +312,6 @@ void recurChangeVisibility(RecurElement *re, const std::vector<const gLevelset *
for(int i = 0; i < re->nbSub(); i++)
recurChangeVisibility(re->sub[i], RPN, TOL);
}
- printf("v=%g vm=%g vis=%d \n",v,vm,re->visible);
}
}
@@ -310,6 +361,7 @@ RecurElement::~RecurElement ()
{
for(int i = 0; i < nbSub(); i++)
if(sub[i]) delete sub[i];
+ delete [] sub;
if(el) delete el;
}
@@ -337,7 +389,7 @@ RecurElement* RecurElement::root()
return this;
}
-double RecurElement::Ls() const
+double RecurElement::ls() const
{
double Tot = 0;
for(int i = 0; i < nbVert(); i++)
@@ -345,292 +397,20 @@ double RecurElement::Ls() const
return Tot / nbVert();
}
-bool RecurElement::cut(int maxlevel, const DI_Element *e, const gLevelset &LS, double TOL,
- std::map<int, double> nodeLs[8])
+bool RecurElement::cut(int maxlevel, const DI_Element *e, std::vector<const gLevelset *> &LsRPN,
+ double TOL, double **nodeLs)
{
- std::vector<const gLevelset *> RPN;
- LS.getRPN(RPN);
recurCut(this, maxlevel, 0);
- bool iC = computeIsCrossed(root(), e, RPN, nodeLs);
+ bool iC = computeIsCrossed(root(), e, LsRPN, nodeLs);
if(TOL < 0.)
recurChangeVisibility(root());
else {
- root()->el->addLs(e, LS);
+ root()->el->addLs(e, LsRPN.back());
recurAddLs(root());
- recurChangeVisibility(root(), RPN, TOL);
+ recurChangeVisibility(root(), LsRPN, TOL);
recurClearLs(root());
}
return iC;
}
-
-/*
-//----------------------------------------------------------------------------------------------
-
-
-
-// create sub-triangles recursively (according to a 4-sub-triangles pattern) until maxlevel
-void recurCut(RecurTriangle *t, int maxlevel, int level)
-{
- if(level++ >= maxlevel) return;
-
- // p3
- // p13 p23
- // p1 p12 p2
- DI_Point p1 = t->tri->pt(0);
- DI_Point p2 = t->tri->pt(1);
- DI_Point p3 = t->tri->pt(2);
- DI_Point p12 ((p1.x()+p2.x())*0.5,(p1.y()+p2.y())*0.5,(p1.z()+p2.z())*0.5);
- DI_Point p13 ((p1.x()+p3.x())*0.5,(p1.y()+p3.y())*0.5,(p1.z()+p3.z())*0.5);
- DI_Point p23 ((p2.x()+p3.x())*0.5,(p2.y()+p3.y())*0.5,(p2.z()+p3.z())*0.5);
- RecurTriangle *tr1 = new RecurTriangle(DI_Triangle(p1,p12,p13));
- recurCut(tr1, maxlevel, level);
- RecurTriangle *tr2 = new RecurTriangle(DI_Triangle(p2,p23,p12));
- recurCut(tr2, maxlevel, level);
- RecurTriangle *tr3 = new RecurTriangle(DI_Triangle(p3,p13,p23));
- recurCut(tr3, maxlevel, level);
- RecurTriangle *tr4 = new RecurTriangle(DI_Triangle(p12,p23,p13));
- recurCut(tr4, maxlevel, level);
- tr1->super = t;
- tr2->super = t;
- tr3->super = t;
- tr4->super = t;
- t->sub[0] = tr1;
- t->sub[1] = tr2;
- t->sub[2] = tr3;
- t->sub[3] = tr4;
-}
-
-// return true if the triangle t->tri is crossed or run along by a primitive levelset in RPN and by the final levelset
-// (the levelset is computed with the values at the nodes of the triangle e)
-bool signChange (RecurTriangle *t, const DI_Triangle *e, const std::vector<const gLevelset *> RPN) {
- bool cS=false;
- std::vector<const gLevelset *> RPNi;
- DI_Triangle tt(*(t->root()->tri));
- std::vector<DI_CuttingPoint> cp;
- for(int l=0;l<(int)RPN.size();l++) {
- const gLevelset *Lsi = RPN[l];
- RPNi.push_back(Lsi);
- if(Lsi->isPrimitive()) {
- tt.computeLs(e,*Lsi);
- for(int i=0;i<(int)cp.size();i++)
- cp[i].addLs(&tt);
- t->tri->addLs(&tt);
- t->tri->getCuttingPoints(&tt,RPNi,cp);
- }
- else {
- for(int p=0;p<(int)cp.size();p++)
- cp[p].chooseLs(Lsi);
- t->tri->chooseLs(Lsi);
- }
- }
- if(cp.size()>1) {
- for(int i=0;i<(int)cp.size();i++)
- if(cp[i].ls()==0) {cS=true; break;}
- }
- t->tri->clearLs();
- return cS;
-}
-
-// Set isCrossed to true if a sub RecurTriangle is crossed.
-// If it has no sub RecurTriangle, set isCrossed to true if the triangle is crossed or run along by the levelset
-//(the levelset is computed with the values at the nodes of the triangle e)
-bool computeIsCrossed (RecurTriangle *t, const DI_Triangle *e, const std::vector<const gLevelset *> RPN) {
- if (t->sub[0]==NULL)
- t->isCrossed = signChange(t,e,RPN);
- else {
- bool iC0 = computeIsCrossed(t->sub[0],e,RPN);
- bool iC1 = computeIsCrossed(t->sub[1],e,RPN);
- bool iC2 = computeIsCrossed(t->sub[2],e,RPN);
- bool iC3 = computeIsCrossed(t->sub[3],e,RPN);
- t->isCrossed = (iC0 || iC1 || iC2 || iC3);
- }
- return t->isCrossed;
-}
-
-// Recursively set the visibility to true if the RecurTriangle is crossed and has no sub recurTriangle
-// or if the RecurTriangle is not crossed and the super RecurTriangle is crossed.
-void changeVisibility(RecurTriangle *t){
- bool superIC = true;
- if(t->super)
- superIC = t->super->isCrossed;
- if((t->isCrossed && t->sub[0]==NULL) || (!t->isCrossed && superIC))
- t->setVisibility(true);
-
- if(t->sub[0]){
- changeVisibility(t->sub[0]);
- changeVisibility(t->sub[1]);
- changeVisibility(t->sub[2]);
- changeVisibility(t->sub[3]);
- }
-}
-
-RecurTriangle* RecurTriangle::root()
-{
- if(super) return super->root();
- return this;
-}
-
-bool RecurTriangle::cut(int maxlevel, const DI_Triangle *e, const std::vector<const gLevelset *> RPN, double TOL)
-{
- recurCut(this, maxlevel, 0);
- bool iC=false;
- if(TOL<0.){
- iC = computeIsCrossed(root(),e,RPN);
- changeVisibility(root());
- }
- else{
- iC = computeIsCrossed(root(),e,RPN);
- changeVisibility(root());
- }
- return iC;
-}
-
-void RecurTriangle::pushTriangles (std::vector<DI_Triangle> &vT)
-{
- if(visible)
- vT.push_back(DI_Triangle(*tri));
- if(sub[0]){
- sub[0]->pushTriangles(vT);
- sub[1]->pushTriangles(vT);
- sub[2]->pushTriangles(vT);
- sub[3]->pushTriangles(vT);
- }
-}
-
-//----------------------------------------------------------------------------------------------
-
-// create sub-quads recursively (according to a 4-sub-quads pattern) until maxlevel
-void recurCut(RecurQuad *q, int maxlevel, int level)
-{
- if(level++ >= maxlevel) return;
-
- // p4 p34 p3
- // p14 p1234 p23
- // p1 p12 p2
- DI_Point p1 = q->quad->pt(0);
- DI_Point p2 = q->quad->pt(1);
- DI_Point p3 = q->quad->pt(2);
- DI_Point p4 = q->quad->pt(3);
- DI_Point p12 ((p1.x()+p2.x())*0.5,(p1.y()+p2.y())*0.5,(p1.z()+p2.z())*0.5);
- DI_Point p23 ((p2.x()+p3.x())*0.5,(p2.y()+p3.y())*0.5,(p2.z()+p3.z())*0.5);
- DI_Point p34 ((p3.x()+p4.x())*0.5,(p3.y()+p4.y())*0.5,(p3.z()+p4.z())*0.5);
- DI_Point p14 ((p1.x()+p4.x())*0.5,(p1.y()+p4.y())*0.5,(p1.z()+p4.z())*0.5);
- DI_Point p1234 ((p12.x()+p34.x())*0.5,(p12.y()+p34.y())*0.5,(p12.z()+p34.z())*0.5);
- RecurQuad *q1 = new RecurQuad(DI_Quad(p1,p12,p1234,p14));
- recurCut(q1, maxlevel, level);
- RecurQuad *q2 = new RecurQuad(DI_Quad(p12,p2,p23,p1234));
- recurCut(q2, maxlevel, level);
- RecurQuad *q3 = new RecurQuad(DI_Quad(p1234,p23,p3,p34));
- recurCut(q3, maxlevel, level);
- RecurQuad *q4 = new RecurQuad(DI_Quad(p14,p1234,p34,p4));
- recurCut(q4, maxlevel, level);
- q1->super = q;
- q2->super = q;
- q3->super = q;
- q4->super = q;
- q->sub[0] = q1;
- q->sub[1] = q2;
- q->sub[2] = q3;
- q->sub[3] = q4;
-}
-
-// return true if the quad q->quad is crossed or run along by a primitive levelset in RPN and by the final levelset
-// (the levelset is computed with the values at the nodes of the quad e)
-bool signChange (RecurQuad *q, const DI_Quad *e, const std::vector<const gLevelset *> RPN) {
- bool cS=false;
- std::vector<const gLevelset *> RPNi;
- DI_Quad qq(*(q->root()->quad));
- std::vector<DI_CuttingPoint> cp;
- for(int l=0;l<(int)RPN.size();l++) {
- const gLevelset *Lsi = RPN[l];
- RPNi.push_back(Lsi);
- if(Lsi->isPrimitive()) {
- qq.computeLs(e,*Lsi);
- for(int i=0;i<(int)cp.size();i++)
- cp[i].addLs(&qq);
- q->quad->addLs(&qq);
- q->quad->getCuttingPoints(&qq,RPNi,cp);
- }
- else {
- for(int p=0;p<(int)cp.size();p++)
- cp[p].chooseLs(Lsi);
- q->quad->chooseLs(Lsi);
- }
- }
- if(cp.size()>1) {
- for(int i=0;i<(int)cp.size();i++)
- if(cp[i].ls()==0) {cS=true; break;}
- }
- q->quad->clearLs();
- return cS;
-}
-
-// Set isCrossed to true if a sub RecurQuad is crossed.
-// If it has no sub RecurQuad, set isCrossed to true if the quad is crossed or run along by the levelset
-//(the levelset is computed with the values at the nodes of the quad e)
-bool computeIsCrossed (RecurQuad *q, const DI_Quad *e, const std::vector<const gLevelset *> RPN) {
- if (q->sub[0]==NULL)
- q->isCrossed = signChange(q,e,RPN);
- else {
- bool iC0 = computeIsCrossed(q->sub[0],e,RPN);
- bool iC1 = computeIsCrossed(q->sub[1],e,RPN);
- bool iC2 = computeIsCrossed(q->sub[2],e,RPN);
- bool iC3 = computeIsCrossed(q->sub[3],e,RPN);
- q->isCrossed = (iC0 || iC1 || iC2 || iC3);
- }
- return q->isCrossed;
-}
-
-// Recursively set the visibility to true if the RecurQuad is crossed and has no sub RecurQuad
-// or if the RecurQuad is not crossed and the parent RecurQuad is crossed.
-void changeVisibility(RecurQuad *q){
- bool superIC = true;
- if(q->super)
- superIC = q->super->isCrossed;
- if((q->isCrossed && q->sub[0]==NULL) || (!q->isCrossed && superIC))
- q->setVisibility(true);
-
- if(q->sub[0]){
- changeVisibility(q->sub[0]);
- changeVisibility(q->sub[1]);
- changeVisibility(q->sub[2]);
- changeVisibility(q->sub[3]);
- }
-}
-
-RecurQuad* RecurQuad::root()
-{
- if(super) return super->root();
- return this;
-}
-
-bool RecurQuad::cut(int maxlevel, const DI_Quad *e, const std::vector<const gLevelset *> RPN, double TOL)
-{
- recurCut(this, maxlevel, 0);
- bool iC=false;
- if(TOL<0.){
- iC = computeIsCrossed(root(),e,RPN);
- changeVisibility(root());
- }
- else{
- iC = computeIsCrossed(root(),e,RPN);
- changeVisibility(root());
- }
- return iC;
-}
-
-void RecurQuad::pushQuads (std::vector<DI_Quad> &vQ)
-{
- if(visible)
- vQ.push_back(DI_Quad(*quad));
- if(sub[0]){
- sub[0]->pushQuads(vQ);
- sub[1]->pushQuads(vQ);
- sub[2]->pushQuads(vQ);
- sub[3]->pushQuads(vQ);
- }
-}
-*/
-
#endif
diff --git a/contrib/DiscreteIntegration/recurCut.h b/contrib/DiscreteIntegration/recurCut.h
index 1144606ca2..132d333695 100644
--- a/contrib/DiscreteIntegration/recurCut.h
+++ b/contrib/DiscreteIntegration/recurCut.h
@@ -20,22 +20,22 @@ class RecurElement
// return the reference element at the root of the tree structure
RecurElement* root();
// return a mean value of the levelset in the element
- inline double Ls() const;
+ inline double ls() const;
// creates RecurElements forming a mesh of the DI_Element e with refined elements close to the zero levelset
// return false if the element is not crossed or run along by the levelset
- bool cut (int maxlevel, const DI_Element *e, const gLevelset &LS, double TOL = -1.,
- std::map<int, double> nodeLs[8] = NULL);
+ bool cut (int maxlevel, const DI_Element *e, std::vector<const gLevelset *> &LsRPN, double TOL = -1.,
+ double **nodeLs = NULL);
};
// push the DI_Elements of the visible RecurElements into v
template <class T>
-static void pushSubElements (RecurElement *re, std::vector<T> &v)
+static void pushSubElements (RecurElement *re, std::vector<T *> &v)
{
if(re->visible)
- v.push_back(T(*((T*)re->el)));
+ v.push_back(new T(*((T*)re->el)));
if(re->sub[0]){
for(int i = 0; i < re->nbSub(); i++)
- pushSubElements(re->sub[i], v);
+ pushSubElements(re->sub[i], v);
}
}
--
GitLab