diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 34c8afd7e2cbd0a3c511622b003534361fc9d076..749adbd68bf5638fa9f158560cd3200ab610b46e 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -165,8 +165,10 @@ void MPolygon::_initVertices()
if(dot(n, ni) < 0.) _parts[i]->revert();
}
- std::vector<MEdge> edg;
- for(unsigned int i = 0; i < _parts.size(); i++) {
+ std::vector<MEdge> edg; // outer edges
+ for(int j = 0; j < _parts[0]->getNumEdges(); j++)
+ edg.push_back(_parts[0]->getEdge(j));
+ for(unsigned int i = 1; i < _parts.size(); i++) {
for(int j = 0; j < 3; j++) {
int k;
for(k = edg.size() - 1; k >= 0; k--)
@@ -183,8 +185,8 @@ void MPolygon::_initVertices()
_edges.push_back(edg[0]);
edg.erase(edg.begin());
while(edg.size()) {
+ MVertex *v = _edges[_edges.size() - 1].getVertex(1);
for(unsigned int i = 0; i < edg.size(); i++) {
- MVertex *v = _edges[_edges.size() - 1].getVertex(1);
if(edg[i].getVertex(0) == v) {
_edges.push_back(edg[i]);
edg.erase(edg.begin() + i);
@@ -195,6 +197,11 @@ void MPolygon::_initVertices()
edg.erase(edg.begin() + i);
break;
}
+ if(i == edg.size() - 1) {
+ _edges.push_back(edg[0]);
+ edg.erase(edg.begin());
+ break;
+ }
}
}
@@ -1118,7 +1125,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
gm->getEntities(gmEntities);
std::vector<const gLevelset *> primitives;
- ls->getPrimitives(primitives);
+ ls->getPrimitivesPO(primitives);
int numVert = gm->indexMeshVertices(true);
int nbLs = (cutElem) ? primitives.size() : 1;
fullMatrix<double> verticesLs(nbLs, numVert + 1);
diff --git a/contrib/DiscreteIntegration/DILevelset.cpp b/contrib/DiscreteIntegration/DILevelset.cpp
index 8ecc23a4a4e0e1f33b433f4a3fc3c81db57ed424..db24a1e6338c7c3d46cc86c65d514806e6321a1a 100644
--- a/contrib/DiscreteIntegration/DILevelset.cpp
+++ b/contrib/DiscreteIntegration/DILevelset.cpp
@@ -49,6 +49,34 @@ void gLevelset::getPrimitives(std::vector<const gLevelset *> &gLsPrimitives) con
}
}
}
+// extrude a list of the primitive levelsets with a "post-order traversal sequence"
+void gLevelset::getPrimitivesPO(std::vector<const gLevelset *> &gLsPrimitives) const {
+ std::stack<const gLevelset *> S;
+ std::stack<const gLevelset *> Sc; // levelset checked
+ S.push(this);
+ while(!S.empty()){
+ const gLevelset *p = S.top();
+ std::vector<const gLevelset *> pp;
+ pp = p->getChildren();
+ if(pp.empty()) {
+ gLsPrimitives.push_back(p);
+ S.pop();
+ }
+ else {
+ if(Sc.empty() || p != Sc.top()) {
+ for(int i = 1; i < (int)pp.size(); i++) Sc.push(p);
+ for(int i = (int)pp.size() - 1; i >= 0; i--) {
+ S.push(pp[i]);
+ if(i > 1) S.push(p);
+ }
+ }
+ else { // levelset has been checked
+ S.pop();
+ Sc.pop();
+ }
+ }
+ }
+}
// return a list with the levelsets in a "Reverse Polish Notation"
void gLevelset::getRPN(std::vector<const gLevelset *> &gLsRPN) const {
@@ -82,7 +110,7 @@ void gLevelset::getRPN(std::vector<const gLevelset *> &gLsRPN) const {
gLevelset::gLevelset(const gLevelset &lv)
{
- tag_=lv.tag_;
+ tag_ = lv.tag_;
}
gLevelsetPlane::gLevelsetPlane(const double * pt, const double *norm, int &tag) : gLevelsetPrimitive(tag) {
@@ -114,13 +142,13 @@ gLevelsetPlane::gLevelsetPlane(const gLevelsetPlane &lv) : gLevelsetPrimitive(lv
x^T A x + [b^T - 2 A t] x + [c - b^T t + t^T A t ] = 0
*/
-gLevelsetQuadric::gLevelsetQuadric(const gLevelsetQuadric &lv):gLevelsetPrimitive(lv)
-{
- for(int i = 0; i < 3; i++){
- B[i] = lv.B[i];
- for(int j = 0; j < 3; j++) A[i][j] = lv.A[i][j];
- }
- C=lv.C;
+gLevelsetQuadric::gLevelsetQuadric(const gLevelsetQuadric &lv) : gLevelsetPrimitive(lv){
+ for(int i = 0; i < 3; i++){
+ B[i] = lv.B[i];
+ for(int j = 0; j < 3; j++)
+ A[i][j] = lv.A[i][j];
+ }
+ C = lv.C;
}
void gLevelsetQuadric::Ax(const double x[3], double res[3], double fact){
for(int i = 0; i < 3; i++){
@@ -133,7 +161,7 @@ void gLevelsetQuadric::Ax(const double x[3], double res[3], double fact){
void gLevelsetQuadric::xAx(const double x[3], double &res, double fact){
res= fact * (A[0][0] * x[0] * x[0] + A[1][1] * x[1] * x[1] + A[2][2] * x[2] * x[2]
- + A[1][0] * x[1] * x[0] * 2. + A[2][0] * x[2] * x[0] * 2. + A[1][2] * x[1] * x[2] * 2.);
+ + A[1][0] * x[1] * x[0] * 2. + A[2][0] * x[2] * x[0] * 2. + A[1][2] * x[1] * x[2] * 2.);
}
void gLevelsetQuadric::translate(const double transl[3]){
@@ -244,7 +272,7 @@ gLevelsetGenCylinder::gLevelsetGenCylinder(const double *pt, const double *dir,
rotate(rot);
translate(pt);
}
-gLevelsetGenCylinder::gLevelsetGenCylinder (const gLevelsetGenCylinder& lv):gLevelsetQuadric(lv){}
+gLevelsetGenCylinder::gLevelsetGenCylinder (const gLevelsetGenCylinder& lv) : gLevelsetQuadric(lv){}
gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, const double &a,
const double &b, const double &c, int &tag) : gLevelsetQuadric(tag) {
@@ -257,7 +285,7 @@ gLevelsetEllipsoid::gLevelsetEllipsoid(const double *pt, const double *dir, cons
rotate(rot);
translate(pt);
}
-gLevelsetEllipsoid::gLevelsetEllipsoid (const gLevelsetEllipsoid& lv):gLevelsetQuadric(lv){}
+gLevelsetEllipsoid::gLevelsetEllipsoid (const gLevelsetEllipsoid& lv) : gLevelsetQuadric(lv){}
gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double &angle, int &tag) : gLevelsetQuadric(tag) {
A[0][0] = 1.;
@@ -269,7 +297,7 @@ gLevelsetCone::gLevelsetCone(const double *pt, const double *dir, const double &
translate(pt);
}
-gLevelsetCone::gLevelsetCone (const gLevelsetCone& lv):gLevelsetQuadric(lv)
+gLevelsetCone::gLevelsetCone (const gLevelsetCone& lv) : gLevelsetQuadric(lv)
{}
gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double *dir, const double &x2, const double &y2, const double &z2,
const double &z, const double &c, int &tag) : gLevelsetQuadric(tag) {
@@ -284,24 +312,23 @@ gLevelsetGeneralQuadric::gLevelsetGeneralQuadric(const double *pt, const double
translate(pt);
}
-gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& lv):gLevelsetQuadric(lv)
+gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric& lv) : gLevelsetQuadric(lv)
{}
-gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv):gLevelset(lv)
+gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv)
{
- std::vector<const gLevelset *> _children=lv.getChildren();
- unsigned siz=_children.size();
- children.resize(siz);
- for(unsigned i=0;i<siz;++i)
- children[i]=_children[i]->clone();
+ std::vector<const gLevelset *> _children=lv.getChildren();
+ unsigned siz = _children.size();
+ children.resize(siz);
+ for(unsigned i = 0; i < siz; ++i)
+ children[i] = _children[i]->clone();
}
-gLevelsetImproved::gLevelsetImproved(const gLevelsetImproved &lv):gLevelset(lv)
-{
- Ls=lv.Ls->clone();
+gLevelsetImproved::gLevelsetImproved(const gLevelsetImproved &lv) : gLevelset(lv){
+ Ls = lv.Ls->clone();
}
gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *dir2, const double *dir3,
- const double &a, const double &b, const double &c, int &tag):gLevelsetImproved() {
+ const double &a, const double &b, const double &c, int &tag) : gLevelsetImproved() {
double dir1m[3] = {-dir1[0], -dir1[1], -dir1[2]};
double dir2m[3] = {-dir2[0], -dir2[1], -dir2[2]};
double dir3m[3] = {-dir3[0], -dir3[1], -dir3[2]};
@@ -321,7 +348,7 @@ gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *d
}
gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *pt3, const double *pt4,
- const double *pt5, const double *pt6, const double *pt7, const double *pt8, int &tag):gLevelsetImproved() {
+ const double *pt5, const double *pt6, const double *pt7, const double *pt8, int &tag) : gLevelsetImproved() {
if(!isPlanar(pt1, pt2, pt3, pt4) || !isPlanar(pt5, pt6, pt7, pt8) || !isPlanar(pt1, pt2, pt5, pt6) ||
!isPlanar(pt3, pt4, pt7, pt8) || !isPlanar(pt1, pt4, pt5, pt8) || !isPlanar(pt2, pt3, pt6, pt7))
printf("WARNING : faces of the box are not planar! %d, %d, %d, %d, %d, %d\n",
@@ -337,9 +364,9 @@ gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *p
Ls = new gLevelsetIntersection(p);
}
-gLevelsetBox::gLevelsetBox(const gLevelsetBox &lv):gLevelsetImproved(lv){}
+gLevelsetBox::gLevelsetBox(const gLevelsetBox &lv) : gLevelsetImproved(lv){}
-gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const double &R, const double &H, int &tag):gLevelsetImproved() {
+gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const double &R, const double &H, int &tag) : gLevelsetImproved() {
double dir2[3] = {-dir[0], -dir[1], -dir[2]};
double n[3]; norm(dir, n);
double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]};
@@ -350,7 +377,7 @@ gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const
Ls = new gLevelsetIntersection(p);
}
-gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const double &R, const double &r, const double &H, int &tag):gLevelsetImproved() {
+gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const double &R, const double &r, const double &H, int &tag) : gLevelsetImproved() {
double dir2[3] = {-dir[0], -dir[1], -dir[2]};
double n[3]; norm(dir, n);
double pt2[3] = {pt[0] + H * n[0], pt[1] + H * n[1], pt[2] + H * n[2]};
@@ -363,12 +390,12 @@ gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const
p2.push_back(new gLevelsetGenCylinder(pt, dir, r, tag));
Ls = new gLevelsetCut(p2);
}
-gLevelsetCylinder::gLevelsetCylinder(const gLevelsetCylinder &lv):gLevelsetImproved(lv){}
+gLevelsetCylinder::gLevelsetCylinder(const gLevelsetCylinder &lv) : gLevelsetImproved(lv){}
gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const double *dir2,
const double &H1, const double &H2, const double &H3,
const double &R1, const double &r1, const double &R2, const double &r2,
- const double &L1, const double &L2, const double &E, int &tag):gLevelsetImproved() {
+ const double &L1, const double &L2, const double &E, int &tag) : gLevelsetImproved() {
double n1[3]; norm(dir1, n1);
double n2[3]; norm(dir2, n2);
double pt1[3] = {pt[0] - n2[0] * H1 / 2., pt[1] - n2[1] * H1 / 2., pt[2] - n2[2] * H1 / 2.};
@@ -398,5 +425,5 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const dou
Ls = new gLevelsetCut(p2);
}
-gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv):gLevelsetImproved(lv){}
+gLevelsetConrod::gLevelsetConrod(const gLevelsetConrod &lv) : gLevelsetImproved(lv){}
#endif
diff --git a/contrib/DiscreteIntegration/DILevelset.h b/contrib/DiscreteIntegration/DILevelset.h
index 1983d76845941e1a897918164224833ce404c83a..848e978eea69055a16b2fa26d5f68dbda64ed1a7 100644
--- a/contrib/DiscreteIntegration/DILevelset.h
+++ b/contrib/DiscreteIntegration/DILevelset.h
@@ -46,6 +46,7 @@ public:
void setTag(int t) { tag_ = t; }
virtual int getTag() const { return tag_; }
void getPrimitives(std::vector<const gLevelset *> &primitives) const;
+ void getPrimitivesPO(std::vector<const gLevelset *> &primitives) const;
void getRPN(std::vector<const gLevelset *> &gLsRPN) const;
double H (const double &x, const double &y, const double &z) const {
if (isInsideDomain(x,y,z) || isOnBorder(x,y,z)) return 1.0;
@@ -143,7 +144,6 @@ public:
gLevelsetQuadric() : gLevelsetPrimitive() {init(); }
gLevelsetQuadric(int &tag) : gLevelsetPrimitive(tag) {init(); }
gLevelsetQuadric(const gLevelsetQuadric &);
-
virtual ~gLevelsetQuadric() {}
double operator () (const double &x, const double &y, const double &z) const;
virtual int type() const = 0;
@@ -226,8 +226,6 @@ public:
if(children.size() != 1) return tag_;
return children[0]->getTag();
}
-
-
};
class gLevelsetReverse : public gLevelset
@@ -254,7 +252,7 @@ public:
gLevelsetCut (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
double choose (double d1, double d2) const {
return (d1 > -d2) ? d1 : -d2; // greater of d1 and -d2
- }
+ }
gLevelsetCut(const gLevelsetCut &lv):gLevelsetTools(lv){}
virtual gLevelset * clone() const{return new gLevelsetCut(*this);}
int type2() const {return CUT;}
@@ -278,10 +276,10 @@ public:
class gLevelsetIntersection : public gLevelsetTools
{
public:
- gLevelsetIntersection (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
- gLevelsetIntersection(const gLevelsetIntersection &lv):gLevelsetTools(lv){}
- virtual gLevelset * clone() const{return new gLevelsetIntersection(*this);}
-
+ gLevelsetIntersection (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+ gLevelsetIntersection(const gLevelsetIntersection &lv):gLevelsetTools(lv) { }
+ virtual gLevelset *clone() const { return new gLevelsetIntersection(*this); }
+
double choose (double d1, double d2) const {
return (d1 > d2) ? d1 : d2; // greater of d1 and d2
}
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 84903e803488dd2f1edc57e18039a749aea0db23..dc7403f22b9db2f6edaae34e77c73ee808f5927a 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -2007,12 +2007,13 @@ bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
pushSubElements(re, ll_subLines);
delete re;
+ int iPrim = 0;
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, l, nodeLs);
+ ll->addLs(this, Lsi, iPrim++, nodeLs);
int nbSubLn = ll_subLines.size();
int nbCP = ll_cp.size();
for(int i = 0; i < nbSubLn; i++) ll_subLines[i]->addLs(ll);
@@ -2039,7 +2040,7 @@ bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- ll->addLs(this, Lsi, l, nodeLs);
+ ll->addLs(this, Lsi, iPrim++, nodeLs);
}
}
}
@@ -2105,12 +2106,13 @@ bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integ
pushSubElements(re, tt_subTriangles);
delete re;
+ int iPrim = 0;
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, l, nodeLs);
+ tt->addLs(this, Lsi, iPrim++, nodeLs);
int nbSubQ = tt_subQuads.size(), nbSubQ1 = nbSubQ;
int nbSubTr = tt_subTriangles.size(), nbSubTr1 = nbSubTr;
int nbSurfLn = tt_surfLines.size(), nbSurfLn1 = nbSurfLn;
@@ -2193,7 +2195,7 @@ bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integ
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- tt->addLs(this, Lsi, l, nodeLs);
+ tt->addLs(this, Lsi, iPrim++, nodeLs);
}
}
}
@@ -2283,12 +2285,13 @@ bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
pushSubElements(re, qq_subQuads);
delete re;
+ int iPrim = 0;
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, l, nodeLs);
+ qq->addLs(this, Lsi, iPrim++, nodeLs);
int nbSubQ = qq_subQuads.size(), nbSubQ1 = nbSubQ;
int nbSubTr = qq_subTriangles.size(), nbSubTr1 = nbSubTr;
int nbSurfLn = qq_surfLines.size(), nbSurfLn1 = nbSurfLn;
@@ -2375,7 +2378,7 @@ bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- qq->addLs(this, Lsi, l, nodeLs);
+ qq->addLs(this, Lsi, iPrim++, nodeLs);
}
}
}
@@ -2479,12 +2482,13 @@ bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrat
pushSubElements(re, tt_subTetras);
delete re;
+ int iPrim = 0;
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, l, nodeLs);
+ tt->addLs(this, Lsi, iPrim++, nodeLs);
int nbSubT = tt_subTetras.size(), nbSubT1 = nbSubT;
int nbSurfQ = tt_surfQuads.size();
int nbSurfTr = tt_surfTriangles.size(), nbSurfTr1 = nbSurfTr;
@@ -2549,7 +2553,7 @@ bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrat
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- tt->addLs(this, Lsi, l, nodeLs);
+ tt->addLs(this, Lsi, iPrim++, nodeLs);
}
}
}
@@ -2652,12 +2656,13 @@ bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
pushSubElements(re, hh_subHexas);
delete re;
+ int iPrim = 0;
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, l, nodeLs);
+ hh->addLs(this, Lsi, iPrim++, nodeLs);
int nbSubH = hh_subHexas.size();
int nbSubT = hh_subTetras.size(), nbSubT1 = nbSubT;
int nbSurfQ = hh_surfQuads.size(), nbSurfQ1 = nbSurfQ;
@@ -2779,7 +2784,7 @@ bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
- hh->addLs(this, Lsi, l, nodeLs);
+ hh->addLs(this, Lsi, iPrim++, nodeLs);
}
}
}