diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 14243e6d6c84bd432fbb18f75cf51fcecefe37d8..d729b27519de6534c979485e02c653215f27e904 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1839,6 +1839,17 @@ GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
Msg::Info("Mesh cutting completed (%g s)", Cpu() - t1);
+ //emi debug
+ // std::vector<GEntity*> entities;
+ // cutGM->getEntities(entities);
+ // for(unsigned int i = 0; i < entities.size(); i++){
+ // std::vector<int> phys = entities[i]->getPhysicalEntities();
+ // for (int j= 0; j < phys.size(); j++){
+ // std::string name = cutGM->getPhysicalName(entities[i]->dim(), phys[j]);
+ // printf("dim =%d elem=%d phys =%s \n", entities[i]->dim(),entities[i]->tag(), name.c_str() );
+ // }
+ // }
+
return cutGM;
}
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 597fc748e1e7048e0e3765748e757c1a3fc20414..ec300ed0ff7646e0a3317430d643a75958ba91d3 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -9,6 +9,7 @@
#include "MElement.h"
#include "MElementCut.h"
#include "gmshLevelset.h"
+#include "MQuadrangle.h"
#if defined(HAVE_DINTEGRATION)
#include "Integration3D.h"
@@ -608,7 +609,7 @@ static int getElementVertexNum(DI_Point *p, MElement *e)
typedef std::set<MVertex*, MVertexLessThanLexicographic> newVerticesContainer;
-static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
+static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
fullMatrix<double> &verticesLs,
GEntity *ge, GModel *GM, int &numEle,
std::map<int, MVertex*> &vertexMap,
@@ -630,7 +631,8 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
std::vector<DI_Hexa *> &hexas)
{
int elementary = ge->tag();
- int eType = e->getTypeForMSH();
+ int eType = e->getType();
+ int recur = e->getPolynomialOrder()-1;
int ePart = e->getPartition();
MElement *eParent = e->getParent();
std::vector<int> gePhysicals = ge->physicals;
@@ -641,7 +643,9 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
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();
MElement *copy = e->copy(vertexMap, newParents, newDomains);
MElement *parent = eParent ? copy->getParent() : copy;
@@ -649,22 +653,23 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
double **nodeLs = new double*[e->getNumPrimaryVertices()];
switch (eType) {
- case MSH_TET_4 :
- case MSH_HEX_8 :
- case MSH_PYR_5 :
- case MSH_PRI_6 :
- case MSH_POLYH_ :
+ case TYPE_TET :
+ case TYPE_HEX :
+ case TYPE_PYR :
+ case TYPE_PRI :
+ case TYPE_POLYH :
{
- if(eType == MSH_TET_4) {
+ if(eType == TYPE_TET) {
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.setPolynomialOrder(recur+1);
for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- tetras, quads, triangles, 0, nodeLs);
+ tetras, quads, triangles, recur, nodeLs);
}
- else if(eType == MSH_HEX_8){
+ else if(eType == TYPE_HEX){
DI_Hexa H(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(),
@@ -673,67 +678,74 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
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.setPolynomialOrder(recur+1);
for(int i = 0; i < 8; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = H.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder, integOrder,
- hexas, tetras, quads, triangles, lines, 0, nodeLs);
+ hexas, tetras, quads, triangles, lines, recur, nodeLs);
}
- else if(eType == MSH_PRI_6){
+ else if(eType == TYPE_PRI){
DI_Tetra T1(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(5)->x(), e->getVertex(5)->y(), e->getVertex(5)->z());
+ T1.setPolynomialOrder(recur+1);
for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
nodeLs[3] = &verticesLs(0, e->getVertex(5)->getIndex());
bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- tetras, quads, triangles, 0, nodeLs);
+ tetras, quads, triangles, recur, 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());
+ T2.setPolynomialOrder(recur+1);
nodeLs[0] = &verticesLs(0, e->getVertex(0)->getIndex());
nodeLs[1] = &verticesLs(0, e->getVertex(4)->getIndex());
nodeLs[2] = &verticesLs(0, e->getVertex(1)->getIndex());
nodeLs[3] = &verticesLs(0, e->getVertex(5)->getIndex());
bool iC2 = T2.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- tetras, quads, triangles, 0, nodeLs);
+ tetras, quads, triangles, recur, 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());
+ T3.setPolynomialOrder(recur+1);
for(int i = 1; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i+2)->getIndex());
bool iC3 = T3.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- tetras, quads, triangles, 0, nodeLs);
+ tetras, quads, triangles, recur, nodeLs);
isCut = iC1 || iC2 || iC3;
}
- else if(eType == MSH_PYR_5){
+ else if(eType == TYPE_PYR){
DI_Tetra T1(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(4)->x(), e->getVertex(4)->y(), e->getVertex(4)->z());
+ T1.setPolynomialOrder(recur+1);
for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
nodeLs[3] = &verticesLs(0, e->getVertex(4)->getIndex());
bool iC1 = T1.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- tetras, quads, triangles, 0, nodeLs);
+ tetras, quads, triangles, recur, 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());
+ T2.setPolynomialOrder(recur+1);
nodeLs[0] = &verticesLs(0, e->getVertex(0)->getIndex());
for(int i = 1; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i+1)->getIndex());
bool iC2 = T2.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- tetras, quads, triangles, 0, nodeLs);
+ tetras, quads, triangles, recur, nodeLs);
isCut = iC1 || iC2;
}
- else if(eType == MSH_POLYH_){
+ else if(eType == TYPE_POLYH){
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());
+ Tet.setPolynomialOrder(recur+1);
for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
bool iC = Tet.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- tetras, quads, triangles, 0, nodeLs);
+ tetras, quads, triangles, recur, nodeLs);
isCut = (isCut || iC);
}
}
@@ -809,18 +821,23 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
if(eParent) {copy->setParent(NULL, false); delete copy;}
}
else { // no cut
- int reg = getElementaryTag(tetras[nbTe]->lsTag(), elementary, newElemTags[3]);
+ int tag ;
+ if (eType == TYPE_HEX)
+ tag = hexas[nbH]->lsTag();
+ else
+ tag = tetras[nbTe]->lsTag();
+ int reg = getElementaryTag(tag, elementary, newElemTags[3]);
std::vector<int> phys;
- getPhysicalTag(tetras[nbTe]->lsTag(), gePhysicals, phys, newPhysTags[3]);
- if(eType == MSH_TET_4)
+ getPhysicalTag(tag, gePhysicals, phys, newPhysTags[3]);
+ if(eType == TYPE_TET)
elements[4][reg].push_back(copy);
- else if(eType == MSH_HEX_8)
+ else if(eType == TYPE_HEX)
elements[5][reg].push_back(copy);
- else if(eType == MSH_PRI_6)
+ else if(eType == TYPE_PRI)
elements[6][reg].push_back(copy);
- else if(eType == MSH_PYR_5)
+ else if(eType == TYPE_PYR)
elements[7][reg].push_back(copy);
- else if(eType == MSH_POLYH_)
+ else if(eType == TYPE_POLYH)
elements[9][reg].push_back(copy);
assignPhysicals(GM, phys, reg, 3, physicals);
}
@@ -866,38 +883,39 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
}
break;
- case MSH_TRI_3 :
- case MSH_TRI_B :
- case MSH_QUA_4 :
- case MSH_POLYG_ :
- case MSH_POLYG_B :
+ case TYPE_TRI :
+ case TYPE_QUA :
+ case TYPE_POLYG :
{
- if((eType == MSH_TRI_3) | (eType == MSH_TRI_B)) {
+ if( eType == TYPE_TRI ) {
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.setPolynomialOrder(recur+1);
for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = T.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- quads, triangles, lines, 0, nodeLs);
+ quads, triangles, lines, recur, nodeLs);
}
- else if(eType == MSH_QUA_4){
+ else if(eType == TYPE_QUA){
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.setPolynomialOrder(recur+1);
for(int i = 0; i < 4; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
isCut = Q.cut(RPN, ipV, ipS, cp, integOrder,integOrder,integOrder,
- quads, triangles, lines, 0, nodeLs);
+ quads, triangles, lines, recur, nodeLs);
}
- else if(eType == MSH_POLYG_ || eType == MSH_POLYG_B){
+ else if(eType == TYPE_POLYG){
for(int i = 0; i < e->getNumChildren(); i++) {
MElement *t = 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());
+ Tri.setPolynomialOrder(recur+1);
for(int i = 0; i < 3; i++) nodeLs[i] = &verticesLs(0, t->getVertex(i)->getIndex());
bool iC = Tri.cut(RPN, ipV, ipS, cp, integOrder, integOrder, integOrder,
- quads, triangles, lines, 0, nodeLs);
+ quads, triangles, lines, recur, nodeLs);
isCut = (isCut || iC);
}
}
@@ -932,6 +950,31 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
else
poly[1].push_back(mt);
}
+ // //if quads
+ // for (unsigned int i = nbQ; i < quads.size(); i++){
+ // MVertex *mv[4] = {NULL, NULL, NULL, NULL};
+ // for(int j = 0; j < 4; j++){
+ // int numV = getElementVertexNum(triangles[i]->pt(j), e);
+ // if(numV == -1) {
+ // MVertex *newv = new MVertex(quads[i]->x(j), quads[i]->y(j), quads[i]->z(j));
+ // std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
+ // mv[j] = *(it.first);
+ // if (!it.second) newv->deleteLast();
+ // }
+ // else {
+ // std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
+ // if(it == vertexMap.end()) {
+ // mv[j] = new MVertex(quads[i]->x(j), quads[i]->y(j),
+ // quads[i]->z(j), 0, numV);
+ // vertexMap[numV] = mv[j];
+ // }
+ // else mv[j] = it->second;
+ // }
+ // }
+ // MQuadrangle *mq = new MQuadrangle(mv[0], mv[1], mv[2], mv[3], 0, 0);
+ // elements[3][reg].push_back(mq);
+ // }
+
bool own = (eParent && !e->ownsParent()) ? false : true;
if(poly[0].size()) {
int n = (e->getParent()) ? e->getNum() : ++numEle;
@@ -984,14 +1027,19 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
if(eParent) {copy->setParent(NULL, false); delete copy;}
}
else { // no cut
- int reg = getElementaryTag(triangles[nbTr]->lsTag(), elementary, newElemTags[2]);
+ int tag;
+ if (eType == TYPE_QUA)
+ tag = quads[nbQ]->lsTag();
+ else
+ tag = triangles[nbTr]->lsTag();
+ int reg = getElementaryTag(tag, elementary, newElemTags[2]);
std::vector<int> phys;
- getPhysicalTag(triangles[nbTr]->lsTag(), gePhysicals, phys, newPhysTags[2]);
- if((eType == MSH_TRI_3) | (eType == MSH_TRI_B))
+ getPhysicalTag(tag, gePhysicals, phys, newPhysTags[2]);
+ if(eType == TYPE_TRI)
elements[2][reg].push_back(copy);
- else if(eType == MSH_QUA_4)
+ else if(eType == TYPE_QUA)
elements[3][reg].push_back(copy);
- else if(eType == MSH_POLYG_ || eType == MSH_POLYG_B)
+ else if(eType == TYPE_POLYG)
elements[8][reg].push_back(copy);
assignPhysicals(GM, phys, reg, 2, physicals);
for(int i = 0; i < 2; i++)
@@ -1040,14 +1088,13 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
}
break;
- case MSH_LIN_2 :
- case MSH_LIN_B :
- case MSH_LIN_C :
+ case TYPE_LIN :
{
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.setPolynomialOrder(recur+1);
for(int i = 0; i < 2; i++) nodeLs[i] = &verticesLs(0, e->getVertex(i)->getIndex());
- isCut = L.cut(RPN, ipV, cp, integOrder, lines, 0, nodeLs);
+ isCut = L.cut(RPN, ipV, cp, integOrder, lines, recur, nodeLs);
if(isCut) {
bool own = (eParent && !e->ownsParent()) ? false : true;
@@ -1099,7 +1146,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
}
}
break;
- case MSH_PNT :
+ case TYPE_PNT :
{
DI_Point P(e->getVertex(0)->x(), e->getVertex(0)->y(), e->getVertex(0)->z());
P.computeLs(RPN.back());
@@ -1138,7 +1185,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::vector<GEntity*> gmEntities;
gm->getEntities(gmEntities);
- std::vector<const gLevelset *> primitives;
+ std::vector<gLevelset *> primitives;
ls->getPrimitivesPO(primitives);
int primS = primitives.size();
int numVert = gm->indexMeshVertices(true);
@@ -1221,7 +1268,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::vector<DI_Quad *> quads;
std::vector<DI_Tetra *> tetras;
std::vector<DI_Hexa *> hexas;
- std::vector<const gLevelset *> RPN;
+ std::vector<gLevelset *> RPN;
ls->getRPN(RPN);
std::vector<int> lsLineRegs;
for(unsigned int i = 0; i < gmEntities.size(); i++) {
diff --git a/Geo/gmshLevelset.cpp b/Geo/gmshLevelset.cpp
index 763b247dcc4e7807466c80c98ddca317e2616443..ef096b97e15eb25051a6df1e4334d01c05d33a35 100644
--- a/Geo/gmshLevelset.cpp
+++ b/Geo/gmshLevelset.cpp
@@ -63,12 +63,12 @@ inline void printNodes(fullMatrix<double> &myNodes, fullMatrix<double> &surf){
}
// extrude a list of the primitive levelsets with a "Level-order traversal sequence"
-void gLevelset::getPrimitives(std::vector<const gLevelset *> &gLsPrimitives) const {
- std::queue<const gLevelset *> Q;
+void gLevelset::getPrimitives(std::vector<gLevelset *> &gLsPrimitives) {
+ std::queue<gLevelset *> Q;
Q.push(this);
while(!Q.empty()){
- const gLevelset *p = Q.front();
- std::vector<const gLevelset *> pp;
+ gLevelset *p = Q.front();
+ std::vector<gLevelset *> pp;
pp = p->getChildren();
if(pp.empty())
gLsPrimitives.push_back(p);
@@ -80,13 +80,13 @@ 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
+void gLevelset::getPrimitivesPO(std::vector<gLevelset *> &gLsPrimitives) {
+ std::stack<gLevelset *> S;
+ std::stack<gLevelset *> Sc; // levelset checked
S.push(this);
while(!S.empty()){
- const gLevelset *p = S.top();
- std::vector<const gLevelset *> pp;
+ gLevelset *p = S.top();
+ std::vector<gLevelset *> pp;
pp = p->getChildren();
if(pp.empty()) {
gLsPrimitives.push_back(p);
@@ -109,13 +109,13 @@ void gLevelset::getPrimitivesPO(std::vector<const gLevelset *> &gLsPrimitives) c
}
// return a list with the levelsets in a "Reverse Polish Notation"
-void gLevelset::getRPN(std::vector<const gLevelset *> &gLsRPN) const {
- std::stack<const gLevelset *> S;
- std::stack<const gLevelset *> Sc; // levelset checked
+void gLevelset::getRPN(std::vector<gLevelset *> &gLsRPN) {
+ std::stack<gLevelset *> S;
+ std::stack<gLevelset *> Sc; // levelset checked
S.push(this);
while(!S.empty()){
- const gLevelset *p = S.top();
- std::vector<const gLevelset *> pp;
+ gLevelset *p = S.top();
+ std::vector<gLevelset *> pp;
pp = p->getChildren();
if(pp.empty()) {
gLsRPN.push_back(p);
@@ -561,7 +561,7 @@ gLevelsetGeneralQuadric::gLevelsetGeneralQuadric (const gLevelsetGeneralQuadric&
gLevelsetTools::gLevelsetTools(const gLevelsetTools &lv) : gLevelset(lv)
{
- std::vector<const gLevelset *> _children=lv.getChildren();
+ std::vector<gLevelset *> _children=lv.getChildren();
unsigned siz = _children.size();
children.resize(siz);
for(unsigned i = 0; i < siz; ++i)
@@ -581,7 +581,7 @@ gLevelsetBox::gLevelsetBox(const double *pt, const double *dir1, const double *d
double n3[3]; norm(dir3, n3);
double pt2[3] = {pt[0] + a * n1[0] + b * n2[0] + c * n3[0], pt[1] + a * n1[1] + b * n2[1] + c * n3[1],
pt[2] + a * n1[2] + b * n2[2] + c * n3[2]};
- std::vector<const gLevelset *> p;
+ std::vector<gLevelset *> p;
p.push_back(new gLevelsetPlane(pt2, dir3, tag++));
p.push_back(new gLevelsetPlane(pt, dir3m, tag++));
p.push_back(new gLevelsetPlane(pt, dir2m, tag++));
@@ -598,7 +598,7 @@ gLevelsetBox::gLevelsetBox(const double *pt1, const double *pt2, const double *p
printf("WARNING : faces of the box are not planar! %d, %d, %d, %d, %d, %d\n",
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));
- std::vector<const gLevelset *> p;
+ std::vector<gLevelset *> p;
p.push_back(new gLevelsetPlane(pt5, pt6, pt8, tag++));
p.push_back(new gLevelsetPlane(pt1, pt4, pt2, tag++));
p.push_back(new gLevelsetPlane(pt1, pt2, pt5, tag++));
@@ -614,7 +614,7 @@ gLevelsetCylinder::gLevelsetCylinder(const double *pt, const double *dir, const
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]};
- std::vector<const gLevelset *> p;
+ std::vector<gLevelset *> p;
p.push_back(new gLevelsetGenCylinder(pt, dir, R, tag++));
p.push_back(new gLevelsetPlane(pt, dir2, tag++));
p.push_back(new gLevelsetPlane(pt2, dir, tag));
@@ -625,11 +625,11 @@ gLevelsetCylinder::gLevelsetCylinder(const double * pt, const double *dir, const
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]};
- std::vector<const gLevelset *> p1;
+ std::vector<gLevelset *> p1;
p1.push_back(new gLevelsetGenCylinder(pt, dir, R, tag++));
p1.push_back(new gLevelsetPlane(pt, dir2, tag++));
p1.push_back(new gLevelsetPlane(pt2, dir, tag++));
- std::vector<const gLevelset *> p2;
+ std::vector<gLevelset *> p2;
p2.push_back(new gLevelsetIntersection(p1));
p2.push_back(new gLevelsetGenCylinder(pt, dir, r, tag));
Ls = new gLevelsetCut(p2);
@@ -658,11 +658,11 @@ gLevelsetConrod::gLevelsetConrod(const double *pt, const double *dir1, const dou
double pt36[3] = {pt35[0] - n3[0] * L2, pt35[1] - n3[1] * L2, pt35[2] - n3[2] * L2};
double pt37[3] = {pt36[0] + n2[0] * H3, pt36[1] + n2[1] * H3, pt36[2] + n2[2] * H3};
double pt38[3] = {pt37[0] + n3[0] * L2, pt37[1] + n3[1] * L2, pt37[2] + n3[2] * L2};
- std::vector<const gLevelset *> p1;
+ std::vector<gLevelset *> p1;
p1.push_back(new gLevelsetBox(pt31, pt32, pt33, pt34, pt35, pt36, pt37, pt38, tag));
p1.push_back(new gLevelsetCylinder(pt1, dir2, R1, H1, tag+6));
p1.push_back(new gLevelsetCylinder(pt2, dir2, R2, H2, tag+9));
- std::vector<const gLevelset *> p2;
+ std::vector<gLevelset *> p2;
p2.push_back(new gLevelsetUnion(p1));
p2.push_back(new gLevelsetGenCylinder(pt1, dir2, r1, tag+12));
p2.push_back(new gLevelsetGenCylinder(pt2, dir2, r2, tag+13));
diff --git a/Geo/gmshLevelset.h b/Geo/gmshLevelset.h
index cfe45bbec5956198ef37194165c19d4427ca148d..1b750089324f79826f90c8766e8b2351d6216b3e 100644
--- a/Geo/gmshLevelset.h
+++ b/Geo/gmshLevelset.h
@@ -55,15 +55,15 @@ public:
bool isInsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain > 0.;}
bool isOutsideDomain (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) * insideDomain < 0.;}
bool isOnBorder (const double &x, const double &y, const double &z) const {return this->operator()(x,y,z) == 0.;}
- virtual std::vector<const gLevelset *> getChildren() const = 0;
+ virtual std::vector<gLevelset *> getChildren() const = 0;
virtual double choose (double d1, double d2) const = 0;
virtual int type() const = 0;
virtual bool isPrimitive() const = 0;
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;
+ void getPrimitives(std::vector<gLevelset *> &primitives) ;
+ void getPrimitivesPO(std::vector<gLevelset *> &primitives) ;
+ void getRPN(std::vector<gLevelset *> &gLsRPN) ;
double H (const double &x, const double &y, const double &z) const {
if (isInsideDomain(x,y,z) || isOnBorder(x,y,z)) return 1.0;
return 0.0;
@@ -105,7 +105,7 @@ public:
tag_ = tag;
}
virtual double operator () (const double x, const double y, const double z) const = 0;
- std::vector<const gLevelset *> getChildren() const { std::vector<const gLevelset *> p; return p; }
+ std::vector<gLevelset *> getChildren() const { std::vector<gLevelset *> p; return p; }
double choose (double d1, double d2) const {
printf("Cannot use function \"choose\" with a primitive!\n");
return d1;
@@ -270,10 +270,10 @@ public:
class gLevelsetTools : public gLevelset
{
protected:
- std::vector<const gLevelset *> children;
+ std::vector<gLevelset *> children;
public:
gLevelsetTools () {}
- gLevelsetTools (std::vector<const gLevelset *> &p) {children = p;}
+ gLevelsetTools (std::vector<gLevelset *> &p) {children = p;}
gLevelsetTools (const gLevelsetTools &);
~gLevelsetTools () {
for(int i = 0; i < (int)children.size(); i++)
@@ -287,7 +287,7 @@ public:
}
return d;
}
- std::vector<const gLevelset *> getChildren() const {
+ std::vector<gLevelset *> getChildren() const {
if(children.size() != 1) return children;
return children[0]->getChildren();
}
@@ -310,13 +310,13 @@ public:
class gLevelsetReverse : public gLevelset
{
protected:
- const gLevelset *ls;
+ gLevelset *ls;
public:
- gLevelsetReverse (const gLevelset *p) : ls(p){}
+ gLevelsetReverse (gLevelset *p) : ls(p){}
double operator () (const double x, const double y, const double z) const {
return -(*ls)(x, y, z);
}
- std::vector<const gLevelset *> getChildren() const {return ls->getChildren();}
+ std::vector<gLevelset *> getChildren() const {return ls->getChildren();}
bool isPrimitive() const {return ls->isPrimitive();}
virtual double choose (double d1, double d2) const {return -ls->choose(d1,d2);}
virtual int type() const {return ls->type();}
@@ -328,7 +328,7 @@ public:
class gLevelsetCut : public gLevelsetTools
{
public:
- gLevelsetCut (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+ gLevelsetCut (std::vector<gLevelset *> p) : gLevelsetTools(p) { }
double choose (double d1, double d2) const {
return (d1 > -d2) ? d1 : -d2; // greater of d1 and -d2
}
@@ -341,7 +341,7 @@ public:
class gLevelsetUnion : public gLevelsetTools
{
public:
- gLevelsetUnion (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+ gLevelsetUnion (std::vector<gLevelset *> p) : gLevelsetTools(p) { }
gLevelsetUnion(const gLevelsetUnion &lv):gLevelsetTools(lv){}
virtual gLevelset * clone() const{return new gLevelsetUnion(*this);}
@@ -355,8 +355,7 @@ public:
class gLevelsetIntersection : public gLevelsetTools
{
public:
- gLevelsetIntersection (std::vector<gLevelset *> &p){ "Coucou here \n";};
- gLevelsetIntersection (std::vector<const gLevelset *> &p) : gLevelsetTools(p) { }
+ gLevelsetIntersection (std::vector<gLevelset *> p) : gLevelsetTools(p) { }
gLevelsetIntersection(const gLevelsetIntersection &lv):gLevelsetTools(lv) { }
virtual gLevelset *clone() const { return new gLevelsetIntersection(*this); }
@@ -370,7 +369,7 @@ public:
class gLevelsetCrack : public gLevelsetTools
{
public:
- gLevelsetCrack (std::vector<const gLevelset *> &p) {
+ gLevelsetCrack (std::vector<gLevelset *> p) {
if (p.size() != 2)
printf("Error : gLevelsetCrack needs 2 levelsets\n");
children.push_back(p[0]);
@@ -393,7 +392,7 @@ public:
gLevelsetImproved(){}
gLevelsetImproved(const gLevelsetImproved &lv);
double operator() (const double x, const double y, const double z) const {return (*Ls)(x, y, z);}
- std::vector<const gLevelset *> getChildren() const { return Ls->getChildren(); }
+ std::vector<gLevelset *> getChildren() const { return Ls->getChildren(); }
double choose (double d1, double d2) const { return Ls->choose(d1, d2); }
virtual int type() const = 0;
bool isPrimitive() const {return Ls->isPrimitive();}
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index df7f99d718975327d0b5607715cc304a191a8f4f..355ead81217ccb57662414d4a7e87e0bfda2dab3 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -6379,7 +6379,7 @@ yyreduce:
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr((yyvsp[(7) - (8)].l)); i++) {
double d; List_Read((yyvsp[(7) - (8)].l), i, &d);
LevelSet *pl = FindLevelSet((int)d);
@@ -6397,7 +6397,7 @@ yyreduce:
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr((yyvsp[(7) - (8)].l)); i++) {
double d; List_Read((yyvsp[(7) - (8)].l), i, &d);
LevelSet *pl = FindLevelSet((int)d);
@@ -6415,7 +6415,7 @@ yyreduce:
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr((yyvsp[(7) - (8)].l)); i++) {
double d; List_Read((yyvsp[(7) - (8)].l), i, &d);
LevelSet *pl = FindLevelSet((int)d);
@@ -6433,7 +6433,7 @@ yyreduce:
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr((yyvsp[(7) - (8)].l)); i++) {
double d; List_Read((yyvsp[(7) - (8)].l), i, &d);
LevelSet *pl = FindLevelSet((int)d);
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 4303ae4228e58599b874fa0224b6123569a04dbe..5f3f1726a1e93748fe9a48f16a63f66388697edc 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1980,7 +1980,7 @@ LevelSet :
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr($7); i++) {
double d; List_Read($7, i, &d);
LevelSet *pl = FindLevelSet((int)d);
@@ -1998,7 +1998,7 @@ LevelSet :
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr($7); i++) {
double d; List_Read($7, i, &d);
LevelSet *pl = FindLevelSet((int)d);
@@ -2016,7 +2016,7 @@ LevelSet :
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr($7); i++) {
double d; List_Read($7, i, &d);
LevelSet *pl = FindLevelSet((int)d);
@@ -2034,7 +2034,7 @@ LevelSet :
yymsg(0, "Levelset %d already exists", t);
}
else {
- std::vector<const gLevelset *> vl;
+ std::vector<gLevelset *> vl;
for(int i = 0; i < List_Nbr($7); i++) {
double d; List_Read($7, i, &d);
LevelSet *pl = FindLevelSet((int)d);
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index c25e70f72ef17cebb8ba7836bde10ac302ff8178..b8d68b5b6443243af97877982ec447193ea4a895 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -64,13 +64,13 @@ 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) {
+ const DI_Element *e, const std::vector<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) {
+ const DI_Element *e, const std::vector<gLevelset *> &RPN) {
double x[3] = {0., 0., 0.};
int n;
for (n = 0; n < el->nbVert(); n++) {
@@ -447,7 +447,7 @@ int bestQuality (DI_Point *p0, DI_Point *p1, DI_Point *p2,
// 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) {
+ const DI_Element *e, const std::vector<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()));
@@ -471,7 +471,7 @@ DI_Point* Newton(const DI_Point *p1, const DI_Point *p2,
// 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) {
+ const DI_Element *e, const std::vector<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.);
midEN.addLs(e);
@@ -522,7 +522,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<gLevelset*> &RPNi) {
int prim = 0; Ls.clear();
if(e->sizeLs() == 0) return;
for(int l = 0; l < (int)RPNi.size(); l++) {
@@ -541,7 +541,7 @@ bool DI_Point::equal(const DI_Point *p) const {
}
// 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<gLevelset*> &RPN) {
int prim = 0; std::vector<double> Ls;
for(int l = 0; l < (int)RPN.size(); l++) {
if(RPN[l]->isPrimitive())
@@ -635,7 +635,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<gLevelset *> &RPNi) {
if(polOrder_ == o) return;
delete [] mid_;
polOrder_ = o;
@@ -720,7 +720,7 @@ 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<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();
@@ -814,7 +814,7 @@ 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,
+ const DI_Element *e, const std::vector<gLevelset *> &RPN,
std::vector<DI_IntegrationPoint *> &ip) const {
std::vector<DI_IntegrationPoint *> ip_ref;
getRefIntegrationPoints(polynomialOrder, ip_ref);
@@ -834,7 +834,7 @@ void DI_Element::integrationPoints (const int polynomialOrder, const DI_Element
ip.push_back(ip_ref[j]);
}
}
-void DI_Element::getCuttingPoints (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+void DI_Element::getCuttingPoints (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_CuttingPoint *> &cp) const {
int s1, s2;
for(int i = 0; i < nbEdg(); i++){
@@ -921,7 +921,7 @@ 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<gLevelset *> &RPN) {
for(int i = 0; i < nbVert(); i++) {
if(pt(i)->isOutsideDomain())
return;
@@ -1412,7 +1412,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<gLevelset *> &RPNi,
std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
{
if (!isCrossed(pt(0), pt(1))) {
@@ -1443,7 +1443,7 @@ void DI_Triangle::splitIntoSubTriangles (std::vector<DI_Triangle *> &subTriangle
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,
+void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cuttingPoints) const
{
@@ -1561,7 +1561,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<gLevelset *> &RPNi,
std::vector<DI_Tetra *> &subTetras, std::vector<DI_Triangle *> &surfTriangles,
std::vector<DI_CuttingPoint *> &cuttingPoints, std::vector<DI_QualError *> &QError) const
{
@@ -1990,12 +1990,12 @@ void DI_Hexa::getRefIntegrationPoints ( const int polynomialOrder , std::vector<
// -----------------------------------------------------------------------------
// cut a line into sublines along the levelset curve
-bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+bool DI_Line::cut (std::vector<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 *> RPNi; RPNi.reserve(RPN.size());
+ std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
DI_Line *ll = new DI_Line(-1, 0, 0, 1, 0, 0, 2.); //reference line
ll->setPolynomialOrder(getPolynomialOrder());
@@ -2010,7 +2010,7 @@ bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
int iPrim = 0;
if(signChange){
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
ll->addLs(this, Lsi, iPrim++, nodeLs);
@@ -2036,7 +2036,7 @@ bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
else{
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
if(Lsi->isPrimitive()) {
ll->addLs(this, Lsi, iPrim++, nodeLs);
}
@@ -2064,7 +2064,7 @@ bool DI_Line::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
// 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,
+bool DI_Line::cut(const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint *> &cp)
{
// check if the line is cut by the level set
@@ -2085,14 +2085,14 @@ bool DI_Line::cut(const DI_Element *e, const std::vector<const gLevelset *> &RPN
}
// cut a triangle into subtriangles along the levelset curve
-bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+bool DI_Triangle::cut (std::vector<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, double **nodeLs) const
{
//printf(" ");print();
- std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
+ std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
DI_Triangle *tt = new DI_Triangle(0, 0, 0, 1, 0, 0, 0, 1, 0, 0.5); //reference triangle
tt->setPolynomialOrder(getPolynomialOrder());
@@ -2109,7 +2109,7 @@ bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integ
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));
+ 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, iPrim++, nodeLs);
@@ -2191,7 +2191,7 @@ bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integ
}
else{
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
if(Lsi->isPrimitive()) {
tt->addLs(this, Lsi, iPrim++, nodeLs);
}
@@ -2236,7 +2236,7 @@ bool DI_Triangle::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integ
}
// cut a triangle into subtriangles along one levelset primitive
-bool DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+bool DI_Triangle::cut (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cp)
{
@@ -2264,14 +2264,14 @@ bool DI_Triangle::cut (const DI_Element *e, const std::vector<const gLevelset *>
}
// cut a quadrangle into subtriangles along the levelset curve
-bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+bool DI_Quad::cut (std::vector<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, double **nodeLs) const
{
//printf(" "); printls();
- std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
+ std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
DI_Quad *qq = new DI_Quad(-1, -1, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, 4.); // reference quad
qq->setPolynomialOrder(getPolynomialOrder());
@@ -2288,7 +2288,7 @@ bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
int iPrim = 0;
if(signChange) {
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
qq->addLs(this, Lsi, iPrim++, nodeLs);
@@ -2374,7 +2374,7 @@ bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
else {
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
if(Lsi->isPrimitive()) {
qq->addLs(this, Lsi, iPrim++, nodeLs);
}
@@ -2419,7 +2419,7 @@ bool DI_Quad::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
// cut a quadrangle into subtriangles along the levelset curve
-bool DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+bool DI_Quad::cut (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_Quad *> &subQuads, std::vector<DI_Triangle *> &subTriangles,
std::vector<DI_Line *> &surfLines, std::vector<DI_CuttingPoint *> &cp)
{
@@ -2458,14 +2458,14 @@ bool DI_Quad::cut (const DI_Element *e, const std::vector<const gLevelset *> &RP
}
// cut a tetrahedron into subtetrahedra along the levelset surface
-bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+bool DI_Tetra::cut (std::vector<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, double **nodeLs) const
{
//printf(" ");print();
- std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
+ std::vector<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());
@@ -2485,7 +2485,7 @@ bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrat
int iPrim = 0;
if(signChange) {
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
tt->addLs(this, Lsi, iPrim++, nodeLs);
@@ -2549,7 +2549,7 @@ bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrat
}
else {
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
if(Lsi->isPrimitive()) {
tt->addLs(this, Lsi, iPrim++, nodeLs);
}
@@ -2600,7 +2600,7 @@ bool DI_Tetra::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrat
}
// cut a tetrahedron into subtetrahedra along the levelset surface
-bool DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+bool DI_Tetra::cut (const DI_Element *e, const std::vector<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)
@@ -2629,7 +2629,7 @@ bool DI_Tetra::cut (const DI_Element *e, const std::vector<const gLevelset *> &R
}
// cut a hex into subtetras along the levelset surface
-bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_IntegrationPoint *> &ip,
+bool DI_Hexa::cut (std::vector<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,
@@ -2638,7 +2638,7 @@ bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
std::vector<DI_Line *> &frontLines, int recurLevel, double **nodeLs) const
{
//printf(" "); print();
- std::vector<const gLevelset *> RPNi; RPNi.reserve(RPN.size());
+ std::vector<gLevelset *> RPNi; RPNi.reserve(RPN.size());
// reference hexa
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.);
@@ -2659,7 +2659,7 @@ bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
int iPrim = 0;
if(signChange){
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
hh->addLs(this, Lsi, iPrim++, nodeLs);
@@ -2780,7 +2780,7 @@ bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
else {
for(int l = 0; l < (int)RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
if(Lsi->isPrimitive()) {
hh->addLs(this, Lsi, iPrim++, nodeLs);
}
@@ -2848,7 +2848,7 @@ bool DI_Hexa::cut (std::vector<const gLevelset *> &RPN, std::vector<DI_Integrati
}
// cut a hex into subtetras along the levelset surface
-bool DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+bool DI_Hexa::cut (const DI_Element *e, const std::vector<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)
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index c2c338fc077050a45928ae62cdd7e19cb87373d7..22d3a6dcdab813b4fdbad4d454c85465a133a479 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -42,7 +42,7 @@ class DI_Point
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<gLevelset *> &RPNi) : x_(x), y_(y), z_(z) {computeLs(e, RPNi);}
// assignment
DI_Point & operator=(const DI_Point & rhs);
// add a levelset value (adjusted to 0 if ls<ZERO_LS_TOL) into the vector Ls
@@ -53,7 +53,7 @@ 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<gLevelset *> &RPNi);
// clear Ls and add the levelset value computed with ls
void computeLs (const gLevelset *ls);
// remove the last value in Ls and add ls
@@ -112,7 +112,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<gLevelset *> &RPNi);
// change the value of weight_
inline void setWeight (double w) {weight_ = w;}
// return the coordinates
@@ -248,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<gLevelset *> &RPNi);
// return tag
int lsTag() const {return lsTag_;}
// return the position of the point with respect to the domain depending on ls_
@@ -279,7 +279,7 @@ 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<gLevelset *> &RPNi);
// return true if the point pt is inside the element
bool contain (const DI_Point *pt) const;
// return true if the element e is inside the element
@@ -303,10 +303,10 @@ class DI_Element
// 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,
+ const std::vector<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,
+ void getCuttingPoints (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_CuttingPoint*> &cp) const;
// return the ith point
virtual DI_Point* pt (int i) const = 0;
@@ -339,7 +339,7 @@ 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<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);
@@ -418,12 +418,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;
- bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ bool cut (std::vector<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,
+ bool cut(const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint*> &cp);
- void selfSplit (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ void selfSplit (const DI_Element *e, const std::vector<gLevelset *> &RPNi,
std::vector<DI_Line *> &subLines, std::vector<DI_CuttingPoint*> &cuttingPoints) const;
inline double length() const {return integral_;}
};
@@ -503,16 +503,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;
- bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ bool cut (std::vector<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, double **nodeLs = NULL) const;
- bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ bool cut (const DI_Element *e, const std::vector<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,
+ void selfSplit (const DI_Element *e, const std::vector<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;
@@ -601,12 +601,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;
- bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ bool cut (std::vector<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, double **nodeLs = NULL) const;
- bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ bool cut (const DI_Element *e, const std::vector<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;
@@ -696,16 +696,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;
- bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ bool cut (std::vector<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, double **nodeLs = NULL) const;
- bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ bool cut (const DI_Element *e, const std::vector<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,
+ void selfSplit ( const DI_Element *e, const std::vector<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;
@@ -825,14 +825,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;
- bool cut (std::vector<const gLevelset *> &LsRPN, std::vector<DI_IntegrationPoint *> &ip,
+ bool cut (std::vector<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, double **nodeLs = NULL) const;
- bool cut (const DI_Element *e, const std::vector<const gLevelset *> &RPNi,
+ bool cut (const DI_Element *e, const std::vector<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);
diff --git a/contrib/DiscreteIntegration/recurCut.cpp b/contrib/DiscreteIntegration/recurCut.cpp
index c9f6320ea5f7cd4c09c490fcdf3462095f7313a8..39a897128c6c2d6f1919587a0b79f3ec1fc1c2a3 100644
--- a/contrib/DiscreteIntegration/recurCut.cpp
+++ b/contrib/DiscreteIntegration/recurCut.cpp
@@ -218,15 +218,15 @@ 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<gLevelset *> &RPN,
double **nodeLs) {
bool cS = false;
DI_Element* elem = re->root()->el;
std::vector<DI_CuttingPoint *> cp;
- std::vector<const gLevelset *> RPNi;
+ std::vector<gLevelset *> RPNi;
int iPrim = 0;
for(unsigned int l = 0; l < RPN.size(); l++) {
- const gLevelset *Lsi = RPN[l];
+ gLevelset *Lsi = RPN[l];
RPNi.push_back(Lsi);
if(Lsi->isPrimitive()) {
elem->addLs(e, Lsi, iPrim++, nodeLs);
@@ -256,7 +256,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<gLevelset *> &RPN,
double **nodeLs) {
if (!re->sub[0])
re->isCrossed = signChange(re, e, RPN, nodeLs);
@@ -284,7 +284,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<gLevelset *> &RPN, double TOL){
printf("rCV : "); re->el->printls();
if(!re->sub[0]){
re->visible = true;
@@ -398,7 +398,7 @@ double RecurElement::ls() const
return Tot / nbVert();
}
-bool RecurElement::cut(int maxlevel, const DI_Element *e, std::vector<const gLevelset *> &LsRPN,
+bool RecurElement::cut(int maxlevel, const DI_Element *e, std::vector<gLevelset *> &LsRPN,
double TOL, double **nodeLs)
{
recurCut(this, maxlevel, 0);
diff --git a/contrib/DiscreteIntegration/recurCut.h b/contrib/DiscreteIntegration/recurCut.h
index aca3b61c4562cdbef859161db71315d208136a5b..569ec66ce90af1777aeba5b6bcd5e64e48319abd 100644
--- a/contrib/DiscreteIntegration/recurCut.h
+++ b/contrib/DiscreteIntegration/recurCut.h
@@ -23,7 +23,7 @@ class RecurElement
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, std::vector<const gLevelset *> &LsRPN, double TOL = -1.,
+ bool cut (int maxlevel, const DI_Element *e, std::vector<gLevelset *> &LsRPN, double TOL = -1.,
double **nodeLs = NULL);
};
diff --git a/gmshpy/gmshGeo.i b/gmshpy/gmshGeo.i
index 8d64807c13d22d3bb23769e55f579c2e284e5746..06c46c7dc093b5e3d82e4ae605a579a8d4e2158f 100644
--- a/gmshpy/gmshGeo.i
+++ b/gmshpy/gmshGeo.i
@@ -49,6 +49,7 @@ namespace std {
%template(GFaceVectorVector) vector< std::vector< GFace *,std::allocator< GFace * > >,std::allocator< std::vector< GFace *,std::allocator< GFace * > > > >;
%template(GFaceList) list<GFace*, std::allocator<GFace*> >;
%template(GEdgeList) list<GEdge*, std::allocator<GEdge*> >;
+ %template(GLevelsetVector) vector<gLevelset *, std::allocator<gLevelset *> >;
}
%include "GmshConfig.h"