diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 44ad5151f0a846db45f22a0e0f190baf446a4e3a..9455547eb79f730feb403eb214a3a89ec900329a 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -461,23 +461,40 @@ static int getBorderTag(int lsTag, int count, int &elementaryMax, std::map<int,
borderTags[lsTag] = reg;
return reg;
}
+ elementaryMax = std::max(elementaryMax, lsTag);
borderTags[lsTag] = lsTag;
return lsTag;
}
-static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
+static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
int &numEle, std::map<int, MVertex*> &vertexMap,
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> newtags[4],
+ 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)
{
int elementary = ge->tag();
int eType = e->getTypeForMSH();
int ePart = e->getPartition();
std::vector<int> gePhysicals = ge->physicals;
+ int integOrder = 1;
+ std::vector<DI_IntegrationPoint> ipV;
+ std::vector<DI_IntegrationPoint> ipS;
+ bool isCut = false;
+ unsigned int nbL = lines.size();
+ unsigned int nbTr = triangles.size();
+ unsigned int nbQ = quads.size();
+ unsigned int nbTe = tetras.size();
+ unsigned int nbH = hexas.size();
+
// copy element
std::vector<MVertex*> vmv;
if(eType != MSH_POLYG_ && eType != MSH_POLYH_) {
@@ -518,23 +535,13 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
case MSH_PRI_6 :
case MSH_POLYH_ :
{
- std::vector<DI_CuttingPoint> cp;
- std::vector<DI_Tetra> subTetras;
- std::vector<DI_Hexa> subHexas;
- std::vector<DI_Quad> surfQuads;
- std::vector<DI_Triangle> surfTriangles;
- std::vector<DI_Line> boundLines;
- int integOrder = 1;
- std::vector<DI_IntegrationPoint> ipV;
- std::vector<DI_IntegrationPoint> ipS;
-
if(eType == MSH_TET_4) {
DI_Tetra 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(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
- T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subTetras, surfQuads, surfTriangles);
+ isCut = T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ tetras, quads, triangles);
}
else if(eType == MSH_HEX_8){
DI_Hexa H(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
@@ -545,8 +552,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
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());
- H.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
- subHexas, subTetras, surfQuads, surfTriangles, boundLines);
+ isCut = H.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
+ hexas, tetras, quads, triangles, lines);
}
else if(eType == MSH_PRI_6){
DI_Tetra T1(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
@@ -561,12 +568,13 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
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());
- T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subTetras, surfQuads, surfTriangles);
- T2.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subTetras, surfQuads, surfTriangles);
- T3.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subTetras, surfQuads, surfTriangles);
+ bool iC1 = T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ tetras, quads, triangles);
+ bool iC2 = T2.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ tetras, quads, triangles);
+ bool iC3 = T3.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ tetras, quads, triangles);
+ isCut = iC1 || iC2 || iC3;
}
else if(eType == MSH_PYR_5){
DI_Tetra T1(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z(),
@@ -577,10 +585,11 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
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());
- T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subTetras, surfQuads, surfTriangles);
- T2.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subTetras, surfQuads, surfTriangles);
+ bool iC1 = T1.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ tetras, quads, triangles);
+ bool iC2 = T2.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ tetras, quads, triangles);
+ isCut = iC1 || iC2;
}
else if(eType == MSH_POLYH_){
for(int i = 0; i < e->getNumChildren(); i++) {
@@ -589,30 +598,27 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
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());
- Tet.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subTetras, surfQuads, surfTriangles);
+ bool iC = Tet.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ tetras, quads, triangles);
+ isCut = (isCut || iC);
}
}
- // if cut
- if((eType == MSH_TET_4 && subTetras.size() > 1) ||
- (eType == MSH_HEX_8 && subTetras.size() > 6) ||
- (eType == MSH_PRI_6 && subTetras.size() > 3) ||
- (eType == MSH_PYR_5 && subTetras.size() > 2) ||
- (eType == MSH_POLYH_ && subTetras.size() > e->getNumChildren())) {
+ MPolyhedron *p1 = NULL, *p2 = NULL;
+ if(isCut) {
std::vector<MTetrahedron*> poly[2];
- for (unsigned int i = 0; i < subTetras.size(); i++){
+ 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(subTetras[i].pt(j), e);
+ int numV = getElementVertexNum(tetras[i].pt(j), e);
if(numV == -1) {
unsigned int k;
for(k = 0; k < newVertices.size(); k++)
- if(equalV(newVertices[k], subTetras[i].pt(j))) break;
+ if(equalV(newVertices[k], tetras[i].pt(j))) break;
if(k == newVertices.size()) {
- mv[j] = new MVertex(subTetras[i].x(j), subTetras[i].y(j),
- subTetras[i].z(j), 0);
+ mv[j] = new MVertex(tetras[i].x(j), tetras[i].y(j),
+ tetras[i].z(j), 0);
newVertices.push_back(mv[j]);
}
else mv[j] = newVertices[k];
@@ -620,65 +626,33 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
if(it == vertexMap.end()) {
- mv[j] = new MVertex(subTetras[i].x(j), subTetras[i].y(j),
- subTetras[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(subTetras[i].lsTag() < 0)
+ if(tetras[i].lsTag() < 0)
poly[0].push_back(mt);
else
poly[1].push_back(mt);
}
- MPolyhedron *p1 = new MPolyhedron(poly[0], copy, true, ++numEle, ePart);
- int reg = getElementaryTag(-1, elementary, newtags[3]);
- elements[9][reg].push_back(p1);
- assignPhysicals(GM, gePhysicals, reg, 3, physicals);
- MPolyhedron *p2 = new MPolyhedron(poly[1], copy, false, ++numEle, ePart);
- elements[9][elementary].push_back(p2);
- assignPhysicals(GM, gePhysicals, elementary, 3, physicals);
-
- for (unsigned int i = 0; i < surfTriangles.size(); i++){
- MVertex *mv[3] = {NULL, NULL, NULL};
- for(int j = 0; j < 3; j++){
- int numV = getElementVertexNum(surfTriangles[i].pt(j), e);
- if(numV == -1) {
- unsigned int k;
- for(k = 0; k < newVertices.size(); k++)
- if(equalV(newVertices[k], surfTriangles[i].pt(j))) break;
- if(k == newVertices.size()) {
- mv[j] = new MVertex(surfTriangles[i].x(j), surfTriangles[i].y(j),
- surfTriangles[i].z(j), 0);
- newVertices.push_back(mv[j]);
- }
- else mv[j] = newVertices[k];
- }
- else {
- std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
- if(it == vertexMap.end()) {
- mv[j] = new MVertex(surfTriangles[i].x(j), surfTriangles[i].y(j),
- surfTriangles[i].z(j), 0, numV);
- vertexMap[numV] = mv[j];
- }
- else mv[j] = it->second;
- }
- }
- MTriangleBorder *tri = new MTriangleBorder(mv[0], mv[1], mv[2],
- p1, p2, ++numEle, ePart);
- double lsT = surfTriangles[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]);
- elements[2][reg].push_back(tri);
- assignPhysicals(GM, gePhysicals, reg, 2, physicals);
+ if(poly[0].size()) {
+ p1 = new MPolyhedron(poly[0], copy, true, ++numEle, ePart);
+ int reg = getElementaryTag(-1, elementary, newtags[3]);
+ elements[9][reg].push_back(p1);
+ assignPhysicals(GM, gePhysicals, reg, 3, physicals);
}
- }
-
+ if(poly[1].size()) {
+ p2 = new MPolyhedron(poly[1], copy, false, ++numEle, ePart);
+ elements[9][elementary].push_back(p2);
+ assignPhysicals(GM, gePhysicals, elementary, 3, physicals);
+ }
+ }
else { // no cut
- int reg = getElementaryTag(subTetras[0].lsTag(), elementary, newtags[3]);
+ int reg = getElementaryTag(tetras[nbTe].lsTag(), elementary, newtags[3]);
if(eType == MSH_TET_4)
elements[4][reg].push_back(copy);
else if(eType == MSH_HEX_8)
@@ -689,38 +663,64 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
elements[7][reg].push_back(copy);
else if(eType == MSH_POLYH_)
elements[9][reg].push_back(copy);
-
assignPhysicals(GM, gePhysicals, reg, 3, physicals);
}
-
+
+ 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);
+ if(numV == -1) {
+ unsigned int k;
+ for(k = 0; k < newVertices.size(); k++)
+ if(equalV(newVertices[k], triangles[i].pt(j))) break;
+ if(k == newVertices.size()) {
+ mv[j] = new MVertex(triangles[i].x(j), triangles[i].y(j),
+ triangles[i].z(j), 0);
+ newVertices.push_back(mv[j]);
+ }
+ else mv[j] = newVertices[k];
+ }
+ 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);
+ vertexMap[numV] = mv[j];
+ }
+ else mv[j] = it->second;
+ }
+ }
+ 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);
+ double 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]);
+ elements[2][reg].push_back(tri);
+ assignPhysicals(GM, gePhysicals, reg, 2, physicals);
+ }
}
break;
case MSH_TRI_3 :
case MSH_QUA_4 :
case MSH_POLYG_ :
{
- std::vector<DI_CuttingPoint> cp;
- std::vector<DI_Quad> subQuads;
- std::vector<DI_Triangle> subTriangles;
- std::vector<DI_Line> boundLines;
- int integOrder = 1;
- std::vector<DI_IntegrationPoint> ipV;
- std::vector<DI_IntegrationPoint> ipS;
-
if(eType == MSH_TRI_3) {
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());
- T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subQuads, subTriangles, boundLines);
+ isCut = T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ quads, triangles, lines);
}
else if(eType == MSH_QUA_4){
DI_Quad Q(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(),
e->getVertex(3)->x(), e->getVertex(3)->y(), e->getVertex(3)->z());
- Q.cut(*ls, ipV, ipS, cp, integOrder,integOrder,integOrder,
- subQuads, subTriangles, boundLines);
+ isCut = Q.cut(*ls, ipV, ipS, cp, integOrder,integOrder,integOrder,
+ quads, triangles, lines);
}
else if(eType == MSH_POLYG_){
for(int i = 0; i < e->getNumChildren(); i++) {
@@ -728,28 +728,27 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
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());
- Tri.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
- subQuads, subTriangles, boundLines);
+ bool iC = Tri.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ quads, triangles, lines);
+ isCut = (isCut || iC);
}
}
- // if cut
- if((eType == MSH_TRI_3 && subTriangles.size() > 1) ||
- (eType == MSH_QUA_4 && subTriangles.size() > 2) ||
- (eType == MSH_POLYG_ && subTriangles.size() > e->getNumChildren())) {
+ MPolygon *p1 = NULL, *p2 = NULL;
+ if(isCut) {
std::vector<MTriangle*> poly[2];
- for (unsigned int i = 0; i < subTriangles.size(); i++){
+ 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(subTriangles[i].pt(j), e);
+ int numV = getElementVertexNum(triangles[i].pt(j), e);
if(numV == -1) {
unsigned int k;
for(k = 0; k < newVertices.size(); k++)
- if(equalV(newVertices[k], subTriangles[i].pt(j))) break;
+ if(equalV(newVertices[k], triangles[i].pt(j))) break;
if(k == newVertices.size()) {
- mv[j] = new MVertex(subTriangles[i].x(j), subTriangles[i].y(j),
- subTriangles[i].z(j), 0);
+ mv[j] = new MVertex(triangles[i].x(j), triangles[i].y(j),
+ triangles[i].z(j), 0);
newVertices.push_back(mv[j]);
}
else mv[j] = newVertices[k];
@@ -757,63 +756,33 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
if(it == vertexMap.end()) {
- mv[j] = new MVertex(subTriangles[i].x(j), subTriangles[i].y(j),
- subTriangles[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(subTriangles[i].lsTag() < 0)
+ if(triangles[i].lsTag() < 0)
poly[0].push_back(mt);
else
poly[1].push_back(mt);
}
- MPolygon *p1 = new MPolygon(poly[0], copy, true, ++numEle, ePart);
- int reg = getElementaryTag(-1, elementary, newtags[2]);
- elements[8][reg].push_back(p1);
- assignPhysicals(GM, gePhysicals, reg, 2, physicals);
- MPolygon *p2 = new MPolygon(poly[1], copy, false, ++numEle, ePart);
- elements[8][elementary].push_back(p2);
- assignPhysicals(GM, gePhysicals, elementary, 2, physicals);
-
- for (unsigned int i = 0; i < boundLines.size(); i++){
- MVertex *mv[2] = {NULL, NULL};
- for(int j = 0; j < 2; j++){
- int numV = getElementVertexNum(boundLines[i].pt(j), e);
- if(numV == -1) {
- unsigned int k;
- for(k = 0; k < newVertices.size(); k++)
- if(equalV(newVertices[k], boundLines[i].pt(j))) break;
- if(k == newVertices.size()) {
- mv[j] = new MVertex(boundLines[i].x(j), boundLines[i].y(j),
- boundLines[i].z(j), 0);
- newVertices.push_back(mv[j]);
- }
- else mv[j] = newVertices[k];
- }
- else {
- std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
- if(it == vertexMap.end()) {
- mv[j] = new MVertex(boundLines[i].x(j), boundLines[i].y(j),
- boundLines[i].z(j), 0, numV);
- vertexMap[numV] = mv[j];
- }
- else mv[j] = it->second;
- }
- }
- MLineBorder *lin = new MLineBorder(mv[0], mv[1], p1, p2, ++numEle, ePart);
- int c = elements[1].count(boundLines[i].lsTag());
- // suppose that the lines have been cut before the surfaces!
- int reg = getBorderTag(boundLines[i].lsTag(), c, newtags[1][0], borderTags[0]);
- elements[1][reg].push_back(lin);
- assignPhysicals(GM, gePhysicals, reg, 1, physicals);
+ if(poly[0].size()) {
+ p1 = new MPolygon(poly[0], copy, true, ++numEle, ePart);
+ int reg = getElementaryTag(-1, elementary, newtags[2]);
+ elements[8][reg].push_back(p1);
+ assignPhysicals(GM, gePhysicals, reg, 2, physicals);
+ }
+ if(poly[1].size()) {
+ p2 = new MPolygon(poly[1], copy, false, ++numEle, ePart);
+ elements[8][elementary].push_back(p2);
+ assignPhysicals(GM, gePhysicals, elementary, 2, physicals);
}
}
-
else { // no cut
- int reg = getElementaryTag(subTriangles[0].lsTag(), elementary, newtags[2]);
+ int reg = getElementaryTag(triangles[nbTr].lsTag(), elementary, newtags[2]);
if(eType == MSH_TRI_3)
elements[2][reg].push_back(copy);
else if(eType == MSH_QUA_4)
@@ -822,20 +791,51 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
elements[8][reg].push_back(copy);
assignPhysicals(GM, gePhysicals, reg, 2, physicals);
}
+
+ 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);
+ if(numV == -1) {
+ unsigned int k;
+ for(k = 0; k < newVertices.size(); k++)
+ if(equalV(newVertices[k], lines[i].pt(j))) break;
+ if(k == newVertices.size()) {
+ mv[j] = new MVertex(lines[i].x(j), lines[i].y(j),
+ lines[i].z(j), 0);
+ newVertices.push_back(mv[j]);
+ }
+ else mv[j] = newVertices[k];
+ }
+ 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);
+ vertexMap[numV] = mv[j];
+ }
+ else mv[j] = it->second;
+ }
+ }
+ 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());
+ // suppose that the lines have been cut before the surfaces!
+ int reg = getBorderTag(lines[i].lsTag(), c, newtags[1][0], borderTags[0]);
+ elements[1][reg].push_back(lin);
+ assignPhysicals(GM, gePhysicals, reg, 1, physicals);
+ }
}
break;
case MSH_LIN_2 :
{
- std::vector<DI_CuttingPoint> cp;
- std::vector<DI_Line> lines;
- int integOrder = 1;
- std::vector<DI_IntegrationPoint> ip;
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());
- L.cut(*ls, ip, cp, integOrder, lines);
+ isCut = L.cut(*ls, ipV, cp, integOrder, lines);
- if(lines.size() > 1) {
- for (unsigned int i = 0; i < lines.size(); i++){
+ 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);
@@ -865,7 +865,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
}
}
else { // no cut
- int reg = getElementaryTag(lines[0].lsTag(), elementary, newtags[1]);
+ int reg = getElementaryTag(lines[nbL].lsTag(), elementary, newtags[1]);
elements[1][reg].push_back(copy);
assignPhysicals(GM, gePhysicals, reg, 1, physicals);
}
@@ -907,15 +907,22 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
for(int i = 0; i < 4; i++)
newtags[i][0] = gm->getMaxElementaryNumber(i);
std::map<int, int> borderTags[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;
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, gmEntities[i], gm, numEle, vertexMap, newVertices,
- elements, physicals, borderTags, newtags);
+ elementCutMesh(e, ls, verticesLs, gmEntities[i], gm, numEle, vertexMap, newVertices,
+ elements, physicals, borderTags, newtags, cp, lines, triangles, quads, tetras, hexas);
cutGM->getMeshPartitions().insert(e->getPartition());
}
+ cp.clear(); lines.clear(); triangles.clear(); quads.clear(); tetras.clear(); hexas.clear();
}
for(unsigned int i = 0; i < newVertices.size(); i++) {
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index bd262edc5653d5433a85513dc85b0565f1cf3b09..c52f486ab741108c2c2fbd03065b9dd31cdb1066 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -263,7 +263,11 @@ class MTriangleBorder : public MTriangle {
}
~MTriangleBorder() {}
MPolyhedron* getDomain(int i) const { return _domains[i]; }
- virtual MElement *getParent() const { return _domains[0]->getParent(); }
+ virtual MElement *getParent() const {
+ if(_domains[0]) return _domains[0]->getParent();
+ if(_domains[1]) return _domains[1]->getParent();
+ return NULL;
+ }
virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
{
return getParent()->getFunctionSpace(order);
@@ -299,7 +303,11 @@ class MLineBorder : public MLine {
}
~MLineBorder() {}
MPolygon* getDomain(int i) const { return _domains[i]; }
- virtual MElement *getParent() const { return _domains[0]->getParent(); }
+ virtual MElement *getParent() const {
+ if(_domains[0]) return _domains[0]->getParent();
+ if(_domains[1]) return _domains[1]->getParent();
+ return NULL;
+ }
virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
{
return getParent()->getFunctionSpace(order);
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index ad888429c3c54586327614ae862bc19b61d69dc7..6219c52786a0ab47c183330a33da6037cf8c8f00 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -1485,7 +1485,7 @@ 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,
+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
{
if (!isCrossed(pt(0), pt(1))) {
@@ -1514,7 +1514,7 @@ void DI_Triangle::splitIntoSubTriangles (std::vector<DI_Triangle> &subTriangles)
subTriangles.push_back(DI_Triangle(pt(2), p02, p12)); // 2-02-12
}
// 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,
+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
{
@@ -1628,7 +1628,7 @@ void DI_Quad::splitIntoTriangles(std::vector<DI_Triangle> &triangles) const
}
// 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,
+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
{
@@ -2043,7 +2043,7 @@ void DI_Hexa::getRefIntegrationPoints ( const int polynomialOrder , std::vector<
// -----------------------------------------------------------------------------
// cut a line into sublines along the levelset curve
-void DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+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) const
{
@@ -2062,7 +2062,7 @@ void DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
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));
+ const gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
ll.addLs(this, *Lsi);
@@ -2110,10 +2110,11 @@ void DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
if(cp[i].equal(ll_cp[p])) {isIn = true; break;}
if(!isIn) cp.push_back(ll_cp[p]);
}
+ 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,
+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
{
// check if the line is cut by the level set
@@ -2133,7 +2134,7 @@ void DI_Line::cut(const DI_Element *e, const std::vector<const gLevelset *> RPNi
}
// cut a triangle into subtriangles along the levelset curve
-void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+bool DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
@@ -2152,7 +2153,7 @@ void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
RecurElement re(&tt);
bool signChange = re.cut(recurLevel, this, Ls);
- re.pushSubTriangles(tt_subTriangles); //for(int i=0;i<(int)tt_subTriangles.size(); i++) tt_subTriangles[i].print();
+ re.pushSubTriangles(tt_subTriangles);
if(signChange){
for(int l = 0; l < (int)RPN.size(); l++) {
@@ -2265,8 +2266,6 @@ void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
}
}*/
- //printf("tt = "); tt.printls();
-
for(int q = 0; q < (int)tt_subQuads.size(); q++) {
tt_subQuads[q].computeLsTagDom(&tt, RPN);
DI_Quad tt_subQ = tt_subQuads[q];
@@ -2282,10 +2281,10 @@ void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
subTriangles.push_back(tt_subTr);
}
for(int l = 0; l < (int)tt_surfLines.size(); l++) {
- tt_surfLines[l].computeLsTagBound(tt_subQuads, tt_subTriangles); //tt_surfLines[l].printls();
+ 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.printls();
+ mappingEl(&tt_surfLn);
tt_surfLn.integrationPoints(polynomialOrderL, &tt_surfLines[l], &tt, RPN, ipS);
surfLines.push_back(tt_surfLn);
}
@@ -2297,10 +2296,11 @@ void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
if(cp[i].equal(tt_cp[p])) {isIn = true; break;}
if(!isIn) cp.push_back(tt_cp[p]);
}
+ 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,
+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
{
@@ -2325,7 +2325,7 @@ void DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *>
}
// cut a quadrangle into subtriangles along the levelset curve
-void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+bool DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
const int polynomialOrderQ, const int polynomialOrderTr, const int polynomialOrderL,
std::vector<DI_Quad> &subQuads, std::vector<DI_Triangle> &subTriangles,
@@ -2346,14 +2346,6 @@ void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
bool signChange = re.cut(recurLevel, this, Ls, 1.e-5);
re.pushSubQuads(qq_subQuads);
- qq.addLs(this, Ls); printf("qq:"); qq.printls();
- for(int i = 0; i < (int)qq_subQuads.size(); i++) {
- qq_subQuads[i].addLs(&qq);
- printf("qi:"); qq_subQuads[i].printls();
- qq_subQuads[i].clearLs();
- }
- qq.clearLs();
-
if(signChange) {
for(int l = 0; l < (int)RPN.size(); l++) {
const gLevelset *Lsi = RPN[l];
@@ -2496,10 +2488,11 @@ void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
if(cp[i].equal(qq_cp[p])) {isIn = true; break;}
if(!isIn) cp.push_back(qq_cp[p]);
}
+ 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,
+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
{
@@ -2535,7 +2528,7 @@ void DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> RPN
}
// cut a tetrahedron into subtetrahedra along the levelset surface
-void DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+bool DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
const int polynomialOrderT, const int polynomialOrderQ, const int polynomialOrderTr,
std::vector<DI_Tetra> &subTetras, std::vector<DI_Quad> &surfQuads,
@@ -2685,10 +2678,11 @@ void DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
if(cp[i].equal(tt_cp[p])) {isIn = true; break;}
if(!isIn) cp.push_back(tt_cp[p]);
}
+ 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,
+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
@@ -2713,7 +2707,7 @@ void DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> RP
}
// cut a hex into subtetras along the levelset surface
-void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+bool DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &cp,
const int polynomialOrderH, const int polynomialOrderT,
const int polynomialOrderQ, const int polynomialOrderTr,
@@ -2937,10 +2931,11 @@ void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
if(cp[i].equal(hh_cp[p])) {isIn = true; break;}
if(!isIn) cp.push_back(hh_cp[p]);
}
+ 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,
+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
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 862e040ab497a96333a1cc1b74638cb9a3b6296d..e31d4ce1812f764e6135ffdf91c8ec8690863133 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -411,12 +411,12 @@ 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;
- void cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+ 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) const;
- void cut(const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ void cut(const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cp) const;
- void selfSplit (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ void selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
std::vector<DI_Line> &subLines, std::vector<DI_CuttingPoint> &cuttingPoints) const;
inline double length() const {return integral_;}
};
@@ -496,16 +496,16 @@ 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;
- void cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+ bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &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) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ 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;
- void selfSplit (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ 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;
double quality () const;
@@ -594,12 +594,12 @@ 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;
- void cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+ bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &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) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ 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;
@@ -689,16 +689,16 @@ 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;
- void cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+ bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &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) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ 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;
- void selfSplit ( const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ 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;
double quality () const;
@@ -818,14 +818,14 @@ 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;
- void cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
+ bool cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
std::vector<DI_IntegrationPoint> &ipS, std::vector<DI_CuttingPoint> &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) const;
- void cut (const DI_Element *e, const std::vector<const gLevelset *> RPNi,
+ 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;
diff --git a/contrib/DiscreteIntegration/recurCut.cpp b/contrib/DiscreteIntegration/recurCut.cpp
index d101fa81879158d483372b041f76fb24bdba38cf..51f24f4b0d3206d07c58cd551b2d5d8fd18994c3 100644
--- a/contrib/DiscreteIntegration/recurCut.cpp
+++ b/contrib/DiscreteIntegration/recurCut.cpp
@@ -166,7 +166,7 @@ void recurCut(RecurElement *re, int maxlevel, int level)
// 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) {
+bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const gLevelset *> &RPN) {
bool cS = false;
DI_Element* elem = re->root()->el;
std::vector<DI_CuttingPoint> cp;
@@ -200,7 +200,7 @@ bool signChange (RecurElement *re, const DI_Element *e, const std::vector<const
// Set isCrossed to true if a sub RecurElement is crossed.
// 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) {
+bool computeIsCrossed (RecurElement *re, const DI_Element *e, const std::vector<const gLevelset *> &RPN) {
if (!re->sub[0])
re->isCrossed = signChange(re, e, RPN);
else {
@@ -227,7 +227,7 @@ void recurChangeVisibility(RecurElement *re){
recurChangeVisibility(re->sub[i]);
}
-void recurChangeVisibility(RecurElement *re, const std::vector<const gLevelset *> RPN, double TOL){
+void recurChangeVisibility(RecurElement *re, const std::vector<const gLevelset *> &RPN, double TOL){
printf("rCV : "); re->el->printls();
if(!re->sub[0]){
re->visible = true;