diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 7e20a17989dd748f41b7b7d5d90f1697953fd4d1..cd129fa7975529f485926e65b20082ec38712b92 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -452,30 +452,43 @@ static void assignPhysicals(GModel *GM, std::vector<int> gePhysicals, int reg, i
}
}
-static int getElementaryTag(double ls, int elementary, std::map<int, int> &newtags)
+static int getElementaryTag(int lsTag, int elementary, std::map<int, int> &newElemTags)
{
- if(ls < 0){
- if(newtags.count(elementary))
- return newtags[elementary];
+ if(lsTag < 0){
+ if(newElemTags.count(elementary))
+ return newElemTags[elementary];
else{
- int reg = ++newtags[0];
- newtags[elementary] = reg;
+ int reg = ++newElemTags[0];
+ newElemTags[elementary] = reg;
return reg;
}
}
return elementary;
}
+static void getPhysicalTag(int lsTag, const std::vector<int> &physicals, std::vector<int> &phys2, std::map<int, int> &newPhysTags)
+{
+ phys2.clear();
+ for(unsigned int i = 0; i < physicals.size(); i++){
+ int phys = physicals[i];
+ if(lsTag < 0){
+ if(!newPhysTags.count(phys))
+ newPhysTags[phys] = ++newPhysTags[0];
+ phys = newPhysTags[phys];
+ }
+ phys2.push_back(phys);
+ }
+}
-static int getBorderTag(int lsTag, int count, int &elementaryMax, std::map<int, int> &borderTags)
+static int getBorderTag(int lsTag, int count, int &maxTag, std::map<int, int> &borderTags)
{
if(borderTags.count(lsTag))
return borderTags[lsTag];
if(count) {
- int reg = ++elementaryMax;
- borderTags[lsTag] = reg;
- return reg;
+ int tag = ++maxTag;
+ borderTags[lsTag] = tag;
+ return tag;
}
- elementaryMax = std::max(elementaryMax, lsTag);
+ maxTag = std::max(maxTag, lsTag);
borderTags[lsTag] = lsTag;
return lsTag;
}
@@ -486,8 +499,10 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
std::vector<MVertex*> &newVertices,
std::map<int, std::vector<MElement*> > elements[10],
std::map<int, std::map<int, std::string> > physicals[4],
- std::map<int, int> borderTags[2],
- std::map<int, int> newtags[4],
+ std::map<int, int> newElemTags[4],
+ 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,
@@ -676,9 +691,11 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
if(poly[0].size()) {
p1 = new MPolyhedron(poly[0], ++numEle, ePart, true, copy);
- int reg = getElementaryTag(-1, elementary, newtags[3]);
+ int reg = getElementaryTag(-1, elementary, newElemTags[3]);
+ std::vector<int> phys;
+ getPhysicalTag(-1, gePhysicals, phys, newPhysTags[3]);
elements[9][reg].push_back(p1);
- assignPhysicals(GM, gePhysicals, reg, 3, physicals);
+ assignPhysicals(GM, phys, reg, 3, physicals);
}
if(poly[1].size()) {
p2 = new MPolyhedron(poly[1], ++numEle, ePart, false, copy);
@@ -687,7 +704,9 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
}
else { // no cut
- int reg = getElementaryTag(tetras[nbTe].lsTag(), elementary, newtags[3]);
+ int reg = getElementaryTag(tetras[nbTe].lsTag(), elementary, newElemTags[3]);
+ std::vector<int> phys;
+ 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)
@@ -698,7 +717,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
elements[7][reg].push_back(copy);
else if(eType == MSH_POLYH_)
elements[9][reg].push_back(copy);
- assignPhysicals(GM, gePhysicals, reg, 3, physicals);
+ assignPhysicals(GM, phys, reg, 3, physicals);
}
for (unsigned int i = nbTr; i < triangles.size(); i++){
@@ -732,9 +751,11 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
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, newtags[2][0], borderTags[1]);
+ int reg = getBorderTag(lsT, c, newElemTags[2][0], borderElemTags[1]);
+ int physTag = getBorderTag(lsT, c, newPhysTags[2][0], borderPhysTags[1]);
+ std::vector<int> phys; phys.push_back(physTag);
elements[2][reg].push_back(tri);
- assignPhysicals(GM, gePhysicals, reg, 2, physicals);
+ assignPhysicals(GM, phys, reg, 2, physicals);
}
}
break;
@@ -812,9 +833,11 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
if(poly[0].size()) {
p1 = new MPolygon(poly[0], ++numEle, ePart, true, copy);
- int reg = getElementaryTag(-1, elementary, newtags[2]);
+ int reg = getElementaryTag(-1, elementary, newElemTags[2]);
+ std::vector<int> phys;
+ getPhysicalTag(-1, gePhysicals, phys, newPhysTags[2]);
elements[8][reg].push_back(p1);
- assignPhysicals(GM, gePhysicals, reg, 2, physicals);
+ assignPhysicals(GM, phys, reg, 2, physicals);
}
if(poly[1].size()) {
p2 = new MPolygon(poly[1], ++numEle, ePart, false, copy);
@@ -823,14 +846,16 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
}
else { // no cut
- int reg = getElementaryTag(triangles[nbTr].lsTag(), elementary, newtags[2]);
+ int reg = getElementaryTag(triangles[nbTr].lsTag(), elementary, newElemTags[2]);
+ std::vector<int> phys;
+ 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)
elements[3][reg].push_back(copy);
else if(eType == MSH_POLYG_)
elements[8][reg].push_back(copy);
- assignPhysicals(GM, gePhysicals, reg, 2, physicals);
+ assignPhysicals(GM, phys, reg, 2, physicals);
}
for (unsigned int i = nbL; i < lines.size(); i++){
@@ -861,11 +886,14 @@ 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 c = elements[1].count(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(lines[i].lsTag(), c, newtags[1][0], borderTags[0]);
+ int reg = getBorderTag(lsL, c, newElemTags[1][0], borderElemTags[0]);
+ int physTag = getBorderTag(lsL, c, newPhysTags[1][0], borderPhysTags[0]);
+ std::vector<int> phys; phys.push_back(physTag);
elements[1][reg].push_back(lin);
- assignPhysicals(GM, gePhysicals, reg, 1, physicals);
+ assignPhysicals(GM, phys, reg, 1, physicals);
}
}
break;
@@ -902,15 +930,19 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
}
}
MLine *ml = new MLine(mv[0], mv[1], ++numEle, ePart);
- int reg = getElementaryTag(lines[i].lsTag(), elementary, newtags[1]);
+ int reg = getElementaryTag(lines[i].lsTag(), elementary, newElemTags[1]);
+ std::vector<int> phys;
+ getPhysicalTag(lines[i].lsTag(), gePhysicals, phys, newPhysTags[1]);
elements[1][reg].push_back(ml);
- assignPhysicals(GM, gePhysicals, reg, 1, physicals);
+ assignPhysicals(GM, phys, reg, 1, physicals);
}
}
else { // no cut
- int reg = getElementaryTag(lines[nbL].lsTag(), elementary, newtags[1]);
+ int reg = getElementaryTag(lines[nbL].lsTag(), elementary, newElemTags[1]);
+ std::vector<int> phys;
+ getPhysicalTag(lines[nbL].lsTag(), gePhysicals, phys, newPhysTags[1]);
elements[1][reg].push_back(copy);
- assignPhysicals(GM, gePhysicals, reg, 1, physicals);
+ assignPhysicals(GM, phys, reg, 1, physicals);
}
}
break;
@@ -918,9 +950,11 @@ static void elementCutMesh(MElement *e, gLevelset *ls,
{
DI_Point P(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z());
P.computeLs(*ls);
- int reg = getElementaryTag(P.lsTag(), elementary, newtags[0]);
+ int reg = getElementaryTag(P.lsTag(), elementary, newElemTags[0]);
+ std::vector<int> phys;
+ getPhysicalTag(P.lsTag(), gePhysicals, phys, newPhysTags[0]);
elements[0][reg].push_back(copy);
- assignPhysicals(GM, gePhysicals, reg, 0, physicals);
+ assignPhysicals(GM, phys, reg, 0, physicals);
}
break;
default :
@@ -946,11 +980,14 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
gm->getEntities(gmEntities);
int numEle = gm->getNumMeshElements();
- std::map<int, int> newtags[4];
- for(int i = 0; i < 4; i++)
- newtags[i][0] = gm->getMaxElementaryNumber(i);
- std::map<int, int> borderTags[2];
-
+ std::map<int, int> newElemTags[4];
+ std::map<int, int> newPhysTags[4];
+ for(int d = 0; d < 4; d++){
+ newElemTags[d][0] = gm->getMaxElementaryNumber(d);
+ newPhysTags[d][0] = gm->getMaxPhysicalNumber(d);
+ }
+ std::map<int, int> borderElemTags[2];
+ std::map<int, int> borderPhysTags[2];
std::vector<const gLevelset *> primitives;
ls->getPrimitives(primitives);
@@ -975,7 +1012,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
MElement *e = gmEntities[i]->getMeshElement(j);
e->setVolumePositive();
elementCutMesh(e, ls, verticesLs, gmEntities[i], gm, numEle, vertexMap, newVertices,
- elements, physicals, borderTags, newtags, cp, lines, triangles, quads, tetras, hexas);
+ elements, physicals, newElemTags, newPhysTags, borderElemTags, borderPhysTags,
+ cp, lines, triangles, quads, tetras, hexas);
cutGM->getMeshPartitions().insert(e->getPartition());
}
cp.clear(); lines.clear(); triangles.clear(); quads.clear(); tetras.clear(); hexas.clear();
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 338dca15e41458404157309a52eb6fb1b6d5b4c8..0ca3b33943cc96d75113f0ea17d9a672d918c9e2 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -189,18 +189,14 @@ int commonV (int &s11, int &s12, int &s21, int &s22) {
return 0;
}
-double adjustLs (double ls) {
- return fabs(ls) < ZERO_LS_TOL ? 0. : ls;
+inline double adjustLs (double ls) {
+ return (fabs(ls) < ZERO_LS_TOL) ? 0. : ls;
}
-bool isCrossed (double ls1, double ls2) {
+inline bool isCrossed (const DI_Point &p1, const DI_Point &p2) {
+ double ls1 = p1.ls(), ls2 = p2.ls();
if(ls1 * ls2 >= 0.) return false;
- if (fabs(ls1) < ZERO_LS_TOL || fabs(ls2) < ZERO_LS_TOL) return false;
return true;
}
-bool isCrossed (const DI_Point &p1, const DI_Point &p2) {
- double v1 = p1.ls(), v2 = p2.ls();
- return isCrossed (v1, v2);
-}
// return the index of the point with minimum x,y and z
int minimum(double *x, double *y, double *z, const int num) {
@@ -523,14 +519,14 @@ void DI_Point::chooseLs (const gLevelset *Lsi) {
double ls1 = Ls[Ls.size() - 2], ls2 = Ls[Ls.size() - 1];
double ls = Lsi->choose(ls1, ls2);
Ls.pop_back(); Ls.pop_back();
- Ls.push_back(ls);
+ addLs(ls);
}
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++) {
if (RPNi[l]->isPrimitive())
- Ls.push_back(e->evalLs(x_, y_, z_, prim++));
+ addLs(e->evalLs(x_, y_, z_, prim++));
else
chooseLs(RPNi[l]);
}
@@ -548,24 +544,35 @@ void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<cons
int prim = 0; std::vector<double> Ls;
for(int l = 0; l < (int)RPN.size(); l++) {
if(RPN[l]->isPrimitive())
- Ls.push_back(e->evalLs(x_, y_, z_, prim++));
+ Ls.push_back(adjustLs(e->evalLs(x_, y_, z_, prim++)));
else {
double ls1 = Ls[Ls.size() - 2], ls2 = Ls[Ls.size() - 1];
double ls = RPN[l]->choose(ls1, ls2);
Ls.pop_back(); Ls.pop_back();
- Ls.push_back(ls);
+ Ls.push_back(adjustLs(ls));
}
}
setLs(Ls.back());
}
// 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_));
+}
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();
- Ls.push_back(ls);
+ 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);
@@ -622,6 +629,9 @@ DI_Element & DI_Element::operator= (const DI_Element &rhs){
return *this;
}
void DI_Element::setPolynomialOrder (int o) {
+ if(polOrder_ == o) return;
+ for(int i = 0; i < nbMid(); i++)
+ delete mid_[i];
polOrder_ = o;
switch (o) {
case 1 :
@@ -642,6 +652,7 @@ void DI_Element::setPolynomialOrder (int o) {
}
}
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];
polOrder_ = o;
@@ -669,10 +680,6 @@ void DI_Element::addLs (const double *ls) {
for(int i = 0; i < nbMid(); i++)
mid_[i]->addLs(ls[nbVert()+i]);
}
-void DI_Element::addLs (int primTag, std::map<int, double> nodeLs[8]) {
- for(int i = 0; i < nbVert(); i++)
- pts_[i]->addLs(nodeLs[i][primTag]);
-}
void DI_Element::addLs (const DI_Element *e) {
if(e->sizeLs() < 1) return;
for(int i = 0; i < nbVert(); i++)
@@ -698,6 +705,24 @@ void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls) {
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)) {
+ addLs(e, Ls);
+ return;
+ }
+ for(int i = 0; i < nbVert(); i++)
+ pts_[i]->addLs(nodeLs[i][primTag]);
+ 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);
+ mid_[j]->addLs(ls);
+ }
+}
void DI_Element::chooseLs (const gLevelset *Lsi) {
if(sizeLs() < 2)
printf("chooseLs with element ls size < 2 : typeEl=%d\n", type());
@@ -708,20 +733,21 @@ void DI_Element::chooseLs (const gLevelset *Lsi) {
}
void DI_Element::clearLs() {
for(int i = 0; i < nbVert(); i++)
- (pts_[i]->Ls).clear();
+ pts_[i]->clearLs();
for(int i = 0; i < nbMid(); i++)
- (mid_[i]->Ls).clear();
+ mid_[i]->clearLs();
}
bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
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;}
+ if(edge >= nbEdg()) {printf("wrong number (%d) for quadratic edge for a ", edge); print(); return false;}
int s1, s2; vert(edge, s1, s2);
- bool quad0 = isQuad();
- if(!quad0) quad(e, RPNi);
+ int order0 = getPolynomialOrder();
+ if(order0 == 1) setPolynomialOrder(2, e, RPNi);
double dist = distance(mid(edge), *xm);
double sideLength = distance(pt(s1), pt(s2));
if (dist / sideLength < 1.e-5) {
- if(!quad0) lin(); printf("dist=%.20f, sideLength=%g, d/sL=%g => do not add quad edge\n", dist, sideLength, dist/sideLength);
+ if(order0 == 1) setPolynomialOrder(1);
+ printf("dist=%.20f, sideLength=%g, d/sL=%g => do not add quad edge\n", dist, sideLength, dist/sideLength);
return true; // do not add the quadratic edge if xm is very close from the middle of the edge
}
DI_Point *p0 = mid_[edge];
@@ -729,12 +755,12 @@ bool DI_Element::addQuadEdge (int edge, DI_Point *xm,
// check if the quadratic edge will produce a twist in the element (detJ<0)
if(!testDetJ()){
// reinitialize quad edges if the element was not quad at the begining
- if(!quad0) lin();
+ if(order0 == 1) setPolynomialOrder(1);
else mid_[edge]->move(p0->x(), p0->y(), p0->z());//mid_[edge] = p0;
printf("detJ<0 when trying to add a quadratic edge in "); print();
return false;
}
- computeIntegral();*/
+ computeIntegral();
return true;
}
bool DI_Element::contain (const DI_Point &p) const {
@@ -863,9 +889,7 @@ double DI_Element::evalLs (const double u, const double v, const double w, int i
}
bool DI_Element::testDetJ() const {
double detJ0 = detJ(x(0), y(0), z(0));
- for(int i = 1; i < nbVert(); i++)
- if(detJ0 * detJ(x(i), y(i), z(i)) < 0.) return false;
- for(int i = nbVert(); i < nbVert() + nbMid(); i++)
+ for(int i = 1; i < nbVert() + nbMid(); i++)
if(detJ0 * detJ(x(i), y(i), z(i)) < 0.) return false;
return true;
}
@@ -889,7 +913,7 @@ double DI_Element::detJ (const double u, const double v, const double w) const {
J[2][0] += x(i) * s[i][2]; J[2][1] += y(i) * s[i][2]; J[2][2] += z(i) * s[i][2];
// delete [] s[i];
}
- delete [] s;
+ delete [] s;
return det3(J[0][0], J[0][1], J[0][2], J[1][0], J[1][1], J[1][2], J[2][0], J[2][1], J[2][2]);
case 2 :
for(int i = 0; i < nbVert() + nbMid(); i++) {
@@ -909,6 +933,7 @@ double DI_Element::detJ (const double u, const double v, const double w) const {
delete [] s;
return sqrt(sq2(J[0][0]) + sq2(J[0][1]) + sq2(J[0][2])); // allways > 0 !!!!
default :
+ delete [] s;
return 1.;
}
}
@@ -1423,12 +1448,12 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevel
std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
std::vector<DI_Line> &surfLines, std::vector<DI_CuttingPoint> &cuttingPoints) const
{
- bool quadT = false; // if true, use quadratic sub triangles
+ 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(fabs(pt(i).ls()) < ZERO_LS_TOL)
+ if(pt(i).isOnBorder())
ze[nbZe++] = i;
for(int i = 0; i < nbZe; i++)
cuttingPoints.push_back(DI_CuttingPoint(pt(ze[i])));
@@ -1538,12 +1563,12 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
std::vector<DI_Tetra> &subTetras, std::vector<DI_Triangle> &surfTriangles,
std::vector<DI_CuttingPoint> &cuttingPoints, std::vector<DI_QualError> &QError) const
{
- bool quadT = false; // use of quadratic surf triangles and quadratic sub tetrahedra
+ 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(fabs(pt(i).ls()) < ZERO_LS_TOL)
+ 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
@@ -1972,8 +1997,7 @@ bool DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) ll.addLs(Lsi->getTag(), nodeLs);
- else ll.addLs(this, *Lsi);
+ ll.addLs(this, *Lsi, nodeLs);
int nbSubLn = ll_subLines.size();
int nbCP = ll_cp.size();
for(int i = 0; i < nbSubLn; i++) ll_subLines[i].addLs(&ll);
@@ -1999,8 +2023,7 @@ bool DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) ll.addLs(Lsi->getTag(), nodeLs);
- else ll.addLs(this, *Lsi);
+ ll.addLs(this, *Lsi, nodeLs);
}
}
}
@@ -2030,7 +2053,7 @@ void DI_Line::cut(const DI_Element *e, const std::vector<const gLevelset *> &RPN
// 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(fabs(pt(i).ls()) < ZERO_LS_TOL) ze[on++] = i;
+ if(pt(i).isOnBorder()) ze[on++] = i;
else if (pt(i).ls() > 0.) pos++;
else neg++;
}
@@ -2070,8 +2093,7 @@ bool DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
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()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
- else tt.addLs(this, *Lsi);
+ tt.addLs(this, *Lsi, nodeLs);
int nbSubQ = tt_subQuads.size(), nbSubQ1 = nbSubQ;
int nbSubTr = tt_subTriangles.size(), nbSubTr1 = nbSubTr;
int nbSurfLn = tt_surfLines.size(), nbSurfLn1 = nbSurfLn;
@@ -2081,7 +2103,7 @@ bool DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
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(fabs(tt.ls(i)) < ZERO_LS_TOL) 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();
@@ -2149,36 +2171,11 @@ bool DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
- else tt.addLs(this, *Lsi);
+ tt.addLs(this, *Lsi, nodeLs);
}
}
}
- //if !changeSign => no cut => clear vectors / triangles.push_back(tt)
- /*if(tt_cp.size() == 1) {
- tt_cp.clear();
- }
- if(tt_cp.size() > 1) {
- bool changeSign = false;
- double sign0 = tt_cp[0].ls();
- for(int i = 1; i < (int)tt_cp.size(); i++)
- if(sign0*tt_cp[i].ls() <= 0) {changeSign = true; break;}
- if(!changeSign) { //print();
- tt_subTriangles.clear();
- tt_subQuads.clear();
- tt_surfLines.clear();
- tt_cp.clear();
- DI_Point vP[3];
- for(int i = 0; i < 3; i++) {
- double ls = Ls(x(i), y(i), z(i));
- vP[i] = DI_Point(tt.x(i), tt.y(i), tt.z(i), ls);
- }
- DI_Triangle trt(vP[0], vP[1], vP[2]); // reference triangle
- tt_subTriangles.push_back(trt);
- }
- }*/
-
for(int q = 0; q < (int)tt_subQuads.size(); q++) {
tt_subQuads[q].computeLsTagDom(&tt, RPN);
DI_Quad tt_subQ = tt_subQuads[q];
@@ -2220,7 +2217,7 @@ void DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *>
// 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(fabs(pt(i).ls()) < ZERO_LS_TOL) ze[on++] = i;
+ if(pt(i).isOnBorder()) ze[on++] = i;
else if (pt(i).ls() > 0.) pos++;
else neg++;
}
@@ -2264,8 +2261,7 @@ bool DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) qq.addLs(Lsi->getTag(), nodeLs);
- else qq.addLs(this, *Lsi);
+ qq.addLs(this, *Lsi, nodeLs);
int nbSubQ = qq_subQuads.size(), nbSubQ1 = nbSubQ;
int nbSubTr = qq_subTriangles.size(), nbSubTr1 = nbSubTr;
int nbSurfLn = qq_surfLines.size(), nbSurfLn1 = nbSurfLn;
@@ -2276,7 +2272,7 @@ bool DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
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(fabs(qq.ls(i)) < ZERO_LS_TOL) ze[cz++] = i;
+ if(qq.pt(i).isOnBorder()) ze[cz++] = i;
else if (qq.ls(i) > 0.) pos++;
else neg++;
}
@@ -2347,33 +2343,11 @@ bool DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) qq.addLs(Lsi->getTag(), nodeLs);
- else qq.addLs(this, *Lsi);
+ qq.addLs(this, *Lsi, nodeLs);
}
}
}
- /*//if !changeSign => clear vectors / quads.push_back(qq)
- if(qq_cp.size() > 1) {
- bool changeSign = false;
- double sign0 = qq_cp[0].ls();
- for(int i = 1; i < (int)qq_cp.size(); i++)
- if(sign0 * qq_cp[i].ls() <= 0) {changeSign = true; break;}
- if(!changeSign) {
- qq_subTriangles.clear();
- qq_subQuads.clear();
- qq_surfLines.clear();
- qq_cp.clear();
- DI_Point vP[4];
- for(int i = 0; i < 4; i++) {
- double ls = Ls(x(i), y(i), z(i));
- vP[i] = DI_Point(qq.x(i), qq.y(i), qq.z(i), ls);
- }
- DI_Quad qt(vP[0], vP[1], vP[2], vP[3]); // reference quad
- qq_subQuads.push_back(qt);
- }
- }*/
-
for(int q = 0; q < (int)qq_subQuads.size(); q++) {
qq_subQuads[q].computeLsTagDom(&qq, RPN);
DI_Quad qq_subQ = qq_subQuads[q];
@@ -2415,7 +2389,7 @@ void DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> &RP
// 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(fabs(pt(i).ls()) < ZERO_LS_TOL) ze[on++] = i;
+ if(pt(i).isOnBorder()) ze[on++] = i;
else if (pt(i).ls() > 0.) pos++;
else neg++;
}
@@ -2473,8 +2447,7 @@ bool DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
- else tt.addLs(this, *Lsi);
+ tt.addLs(this, *Lsi, nodeLs);
int nbSubT = tt_subTetras.size(), nbSubT1 = nbSubT;
int nbSurfQ = tt_surfQuads.size();
int nbSurfTr = tt_surfTriangles.size(), nbSurfTr1 = nbSurfTr;
@@ -2484,7 +2457,7 @@ bool DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
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(fabs(tt.ls(i)) < ZERO_LS_TOL) 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();
@@ -2535,34 +2508,11 @@ bool DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) tt.addLs(Lsi->getTag(), nodeLs);
- else tt.addLs(this, *Lsi);
+ tt.addLs(this, *Lsi, nodeLs);
}
}
}
- /*//if !changeSign => clear vectors / tetras.push_back(tt)
- if(tt_cp.size() > 1) {
- bool changeSign = false;
- double sign0 = tt_cp[0].ls();
- for(int i = 1; i < (int)tt_cp.size(); i++)
- if(sign0 * tt_cp[i].ls() <= 0) {changeSign = true; break;}
- if(!changeSign) {
- tt_subTetras.clear();
- tt_surfQuads.clear();
- tt_surfTriangles.clear();
- tt_cp.clear();
- QError.clear();
- DI_Point vP[4];
- for(int i = 0; i < 4; i++) {
- double ls = Ls(x(i), y(i), z(i));
- vP[i] = DI_Point(tt.x(i), tt.y(i), tt.z(i), ls);
- }
- DI_Tetra tet(vP[0], vP[1], vP[2], vP[3]); // reference tetra
- tt_subTetras.push_back(tet);
- }
- }*/
-
for(int i = 0; i < (int)QError.size(); i++)
QError[i].print(this);
@@ -2609,7 +2559,7 @@ void DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> &R
// 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(fabs(pt(i).ls()) < ZERO_LS_TOL) ze[on++] = i;
+ if(pt(i).isOnBorder()) ze[on++] = i;
else if (pt(i).ls() > 0.) pos++;
else neg++;
}
@@ -2658,8 +2608,7 @@ bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) hh.addLs(Lsi->getTag(), nodeLs);
- else hh.addLs(this, *Lsi);
+ hh.addLs(this, *Lsi, nodeLs);
int nbSubH = hh_subHexas.size();
int nbSubT = hh_subTetras.size(), nbSubT1 = nbSubT;
int nbSurfQ = hh_surfQuads.size(), nbSurfQ1 = nbSurfQ;
@@ -2674,7 +2623,7 @@ bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
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(fabs(hh.ls(i)) < ZERO_LS_TOL) ze[cz++] = i;
+ if(hh.pt(i).isOnBorder()) ze[cz++] = i;
else if (hh.ls(i) > 0.) pos++;
else neg++;
}
@@ -2774,37 +2723,11 @@ bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) hh.addLs(Lsi->getTag(), nodeLs);
- else hh.addLs(this, *Lsi);
+ hh.addLs(this, *Lsi, nodeLs);
}
}
}
- /*//if the sign does not change and !=0 in the hexa => returns only a hexa
- if(hh_cp.size() > 1) {
- bool changeSign = false;
- double sign0 = hh_cp[0].ls();
- for(int i = 1; i < (int)hh_cp.size(); i++)
- if(sign0*hh_cp[i].ls() <= 0)
- {changeSign = true; break;}
- if(!changeSign) {
- hh_subHexas.clear();
- hh_subTetras.clear();
- hh_surfQuads.clear();
- hh_surfTriangles.clear();
- hh_frontLines.clear();
- hh_cp.clear();
- QError.clear();
- DI_Point vP[8];
- for(int i = 0; i < 8; i++) {
- double ls = Ls(x(i), y(i), z(i));
- vP[i] = DI_Point(hh.x(i), hh.y(i), hh.z(i), ls);
- }
- DI_Hexa hex(vP[0], vP[1], vP[2], vP[3], vP[4], vP[5], vP[6], vP[7]); // reference hexa
- hh_subHexas.push_back(hex);
- }
- }*/
-
for(int i = 0; i < (int)QError.size(); i++)
QError[i].print(this);
@@ -2866,7 +2789,7 @@ void DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> &RP
// 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(fabs(pt(i).ls()) < ZERO_LS_TOL) ze[on++] = i;
+ if(pt(i).isOnBorder()) ze[on++] = i;
else if (pt(i).ls() > 0.) pos++;
else neg++;
}
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 8cac2681b26bf35a95b2feeaf700a4d2c8084c08..791fbea19ac7d372faf911da624ce87ed4da497d 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -30,14 +30,14 @@ class DI_Point
{
// coordinates of the point
double x_, y_, z_;
- public :
// vector containing the levelset values of the point
std::vector<double> Ls;
+ public :
// 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, const double ls) : x_(x), y_(y), z_(z) {Ls.push_back(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);}
DI_Point(const DI_Point &p) : x_(p.x()), y_(p.y()), z_(p.z()) {Ls.clear(); Ls = p.Ls;}
@@ -46,9 +46,9 @@ class DI_Point
// destructor
~DI_Point () {Ls.clear();}
// add a levelset value (adjusted to 0 if ls<ZERO_LS_TOL) into the vector Ls
- inline void addLs (const double ls);
+ void addLs (const double ls);
// add a levelset value evaluated into the element e
- inline void addLs (const DI_Element *e);
+ void addLs (const DI_Element *e);
// choose the value of the levelset among the last two levelset values of Ls,
// delete the last two values and add the chosen one
void chooseLs (const gLevelset *Lsi);
@@ -57,7 +57,7 @@ class DI_Point
// clear Ls and add the levelset value computed with ls
void computeLs (const gLevelset &ls);
// remove the last value in Ls and add ls
- inline void changeLs (const double ls) {Ls.pop_back(); Ls.push_back(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)
@@ -70,6 +70,10 @@ class DI_Point
inline double ls() const {return Ls.back();}
// return the ith value of Ls
inline double ls(int i) const {return Ls[i];}
+ // return the size of the vector of ls values
+ inline int sizeLs() const {return Ls.size();}
+ // clear the values in Ls
+ inline void clearLs() {Ls.clear();}
// 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.;}
@@ -80,7 +84,7 @@ class DI_Point
return 0;
}
// print the coordiantes
- void print() const {printf("Point (%g,%g,%g)\n", x_, y_, z_);}
+ inline void print() const {printf("Point (%g,%g,%g)\n", x_, y_, z_);}
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]);
@@ -229,13 +233,13 @@ class DI_Element
virtual double refIntegral() const = 0;
// add a levelset value to each point
void addLs (const double *ls);
- // add the levelset value to each point from the map with each primitive value
- void addLs (int primTag, std::map<int, double> nodeLs[8]);
// evaluate the levelset value at each point and add it to each point
void addLs (const DI_Element *e);
- // add the level set at each vertex of the real element e (same type as this)
+ // 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[]);
// clear the levelset of the vertices
void clearLs();
// compute the levelset values at the mid edge points and add a quadratic edge,
@@ -286,7 +290,7 @@ class DI_Element
inline double ls (int i, int j) const {
return i < nbVert() ? pts_[i]->ls(j) : mid_[i - nbVert()]->ls(j);}
// return the number of levelset values of the points
- inline int sizeLs() const {return pts_[0]->Ls.size();}
+ inline int sizeLs() const {return pts_[0]->sizeLs();}
// return the interpolating nodal shape functions evaluated at point (u,v,w)
// in parametric coordinates (if order = -1, use the polynomial order of the element)
virtual void getShapeFunctions (double u, double v, double w, double s[], int order = -1) const = 0;
@@ -322,12 +326,11 @@ class DI_CuttingPoint
double xl_, yl_, zl_;
std::vector<double> Ls;
public:
- DI_CuttingPoint (const DI_Point &pt)
- : x_(pt.x()), y_(pt.y()), z_(pt.z()), xl_(pt.x()), yl_(pt.y()), zl_(pt.z()), Ls(pt.Ls) { }
+ 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;}
- inline void addLs (const DI_Element *e) { Ls.push_back(e->evalLs(x_, y_, z_));}
- inline void addLs (double ls) {Ls.push_back(ls);}
+ 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_;}
@@ -337,13 +340,13 @@ class DI_CuttingPoint
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
- bool isInsideDomain () const {return Ls.back() < 0.;}
- bool isOutsideDomain () const {return Ls.back() > 0.;}
- bool isOnBorder () const {return Ls.back() == 0.;}
+ 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;
- void print() 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 {
diff --git a/contrib/DiscreteIntegration/recurCut.cpp b/contrib/DiscreteIntegration/recurCut.cpp
index e93eb2d90ed45869e75ebd4bad82fae3d88e9f98..a54f06376d7a3a82159fc4bcb8a38514dd134448 100644
--- a/contrib/DiscreteIntegration/recurCut.cpp
+++ b/contrib/DiscreteIntegration/recurCut.cpp
@@ -176,8 +176,7 @@ bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- if(nodeLs && nodeLs[0].count(Lsi->getTag())) elem->addLs(Lsi->getTag(), nodeLs);
- else elem->addLs(e, *Lsi);
+ elem->addLs(e, *Lsi, nodeLs);
for(unsigned int i = 0; i < cp.size(); i++)
cp[i].addLs(elem);
if (re->super) re->el->addLs(elem);