diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index df42a1b8f59d9a05d52f14b0cd869ac3133e0236..79bdfe1e2225b4252a74f0e8a53f21f3f5a4f6a8 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -319,7 +319,7 @@ int GModel::readMSH(const std::string &name)
fscanf(fp, "%d", &numVertices);
}
}
- int indices[60];
+ int *indices = new int[numVertices];
for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &indices[j]);
std::vector<MVertex*> vertices;
if(vertexVector.size()){
@@ -332,6 +332,7 @@ int GModel::readMSH(const std::string &name)
vertices, elements, physicals);
if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, "Reading elements");
+ delete [] indices;
}
}
else{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 4f90fc0f76e76fe75dbf714e070d8d09b166d2e9..2c70f4cb53641889a86c41ae4b225cd8c8379a4b 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -82,6 +82,10 @@ class MElement
// get the vertices
virtual int getNumVertices() const = 0;
virtual MVertex *getVertex(int num) = 0;
+ virtual void getCoordinates(int num, double c[3]) const
+ {
+ c[0] = 0.; c[1] = 0.; c[2] = 0.;
+ }
// get the vertex using the I-deas UNV ordering
virtual MVertex *getVertexUNV(int num){ return getVertex(num); }
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 4c54c68ac22e09cf95abe5fdcb2087585e06fe84..38cda7e2b039e53c93f1448ec08a4251ebc571b1 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -161,7 +161,7 @@ void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
if(physical < 0) revert();
fprintf(fp, " %d", n);
- int verts[180];
+ int *verts = new int[n];
for(unsigned int i = 0; i < _parts.size(); i++)
for(int j = 0; j < 4; j++)
verts[i * 4 + j] = _parts[i]->getVertex(j)->getIndex();
@@ -176,6 +176,7 @@ void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
}
if(physical < 0) revert();
+ delete [] verts;
}
//------------------------------------------- MPolygon ------------------------------
@@ -332,7 +333,7 @@ void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
if(physical < 0) revert();
fprintf(fp, " %d", n);
- int verts[120];
+ int *verts = new int[n];
for(unsigned int i = 0; i < _parts.size(); i++)
for(int j = 0; j < 3; j++)
verts[i * 3 + j] = _parts[i]->getVertex(j)->getIndex();
@@ -347,6 +348,64 @@ void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
}
if(physical < 0) revert();
+ delete [] verts;
+}
+
+//----------------------------------- MTriangleBorder ------------------------------
+
+void MTriangleBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+{
+ double uvw[3][3];
+ for(int j = 0; j < 3; j++) {
+ double xyz[3]; getCoordinates(j, xyz);
+ getParent()->xyz2uvw(xyz, uvw[j]);
+ }
+ MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
+ MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
+ MVertex v2(uvw[2][0], uvw[2][1], uvw[2][2]);
+ MTriangle tt(&v0, &v1, &v2);
+ tt.getIntegrationPoints(pOrder, npts, pts);
+ double jac[3][3];
+ for(int ip = 0; ip < (*npts); ip++){
+ const double u = pts[ip]->pt[0];
+ const double v = pts[ip]->pt[1];
+ const double w = pts[ip]->pt[2];
+ const double weight = pts[ip]->weight;
+ const double detJ = tt.getJacobian(u, v, w, jac);
+ SPoint3 p; tt.pnt(u, v, w, p);
+ (*pts[ip]).pt[0] = p.x();
+ (*pts[ip]).pt[1] = p.y();
+ (*pts[ip]).pt[2] = p.z();
+ (*pts[ip]).weight = detJ * weight;
+ }
+}
+
+//-------------------------------------- MLineBorder ------------------------------
+
+void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
+{
+ double uvw[2][3];
+ for(int j = 0; j < 2; j++) {
+ double xyz[3]; getCoordinates(j, xyz);
+ getParent()->xyz2uvw(xyz, uvw[j]);
+ }
+ MVertex v0(uvw[0][0], uvw[0][1], uvw[0][2]);
+ MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
+ MLine ll(&v0, &v1);
+ ll.getIntegrationPoints(pOrder, npts, pts);
+ double jac[3][3];
+ for(int ip = 0; ip < (*npts); ip++){
+ const double u = pts[ip]->pt[0];
+ const double v = pts[ip]->pt[1];
+ const double w = pts[ip]->pt[2];
+ const double weight = pts[ip]->weight;
+ const double detJ = ll.getJacobian(u, v, w, jac);
+ SPoint3 p; ll.pnt(u, v, w, p);
+ (*pts[ip]).pt[0] = p.x();
+ (*pts[ip]).pt[1] = p.y();
+ (*pts[ip]).pt[2] = p.z();
+ (*pts[ip]).weight = detJ * weight;
+ }
}
//---------------------------------------- CutMesh ----------------------------
@@ -392,13 +451,25 @@ static int getElementaryTag(double ls, int elementary, std::map<int, int> &newta
return elementary;
}
+static int getBorderTag(int lsTag, int count, int &elementaryMax, std::map<int, int> &borderTags)
+{
+ if(borderTags.count(lsTag))
+ return borderTags[lsTag];
+ if(count) {
+ int reg = ++elementaryMax;
+ borderTags[lsTag] = reg;
+ return reg;
+ }
+ borderTags[lsTag] = lsTag;
+ return lsTag;
+}
+
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::vector<MElement*> > border[2],
std::map<int, std::map<int, std::string> > physicals[4],
- std::map<int, std::vector<GEntity*> > &entityCut,
+ std::map<int, int> borderTags[2],
std::map<int, int> newtags[4])
{
int elementary = ge->tag();
@@ -408,15 +479,32 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
// copy element
std::vector<MVertex*> vmv;
- for(int i = 0; i < e->getNumVertices(); i++) {
- int numV = e->getVertex(i)->getNum();
- if(vertexMap.count(numV))
- vmv.push_back(vertexMap[numV]);
- else {
- MVertex *mv = new MVertex(e->getVertex(i)->x(), e->getVertex(i)->y(),
- e->getVertex(i)->z(), 0, numV);
- vmv.push_back(mv);
- vertexMap[numV] = mv;
+ if(eType != MSH_POLYG_ && eType != MSH_POLYH_) {
+ for(int i = 0; i < e->getNumVertices(); i++) {
+ MVertex *vert = e->getVertex(i);
+ int numV = vert->getNum();
+ if(vertexMap.count(numV))
+ vmv.push_back(vertexMap[numV]);
+ else {
+ MVertex *mv = new MVertex(vert->x(), vert->y(), vert->z(), 0, numV);
+ vmv.push_back(mv);
+ vertexMap[numV] = mv;
+ }
+ }
+ }
+ else {
+ for(int i = 0; i < e->getNumChildren(); i++) {
+ for(int j = 0; j < e->getChild(i)->getNumVertices(); j++) {
+ MVertex *vert = e->getChild(i)->getVertex(j);
+ int numV = vert->getNum();
+ if(vertexMap.count(numV))
+ vmv.push_back(vertexMap[numV]);
+ else {
+ MVertex *mv = new MVertex(vert->x(), vert->y(), vert->z(), 0, numV);
+ vmv.push_back(mv);
+ vertexMap[numV] = mv;
+ }
+ }
}
}
MElementFactory factory;
@@ -427,6 +515,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
case MSH_HEX_8 :
case MSH_PYR_5 :
case MSH_PRI_6 :
+ case MSH_POLYH_ :
{
std::vector<DI_CuttingPoint> cp;
std::vector<DI_Tetra> subTetras;
@@ -492,12 +581,24 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
T2.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
subTetras, surfQuads, surfTriangles);
}
+ else if(eType == MSH_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.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ subTetras, surfQuads, surfTriangles);
+ }
+ }
// 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_PYR_5 && subTetras.size() > 2) ||
+ (eType == MSH_POLYH_ && subTetras.size() > e->getNumChildren())) {
std::vector<MTetrahedron*> poly[2];
for (unsigned int i = 0; i < subTetras.size(); i++){
@@ -566,8 +667,12 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
}
MTriangleBorder *tri = new MTriangleBorder(mv[0], mv[1], mv[2],
p1, p2, ++numEle, ePart);
- border[1][surfTriangles[i].lsTag()].push_back(tri);
- entityCut[surfTriangles[i].lsTag()].push_back(ge);
+ 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);
}
}
@@ -581,6 +686,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
elements[6][reg].push_back(copy);
else if(eType == MSH_PYR_5)
elements[7][reg].push_back(copy);
+ else if(eType == MSH_POLYH_)
+ elements[9][reg].push_back(copy);
assignPhysicals(GM, gePhysicals, reg, 3, physicals);
}
@@ -589,6 +696,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
break;
case MSH_TRI_3 :
case MSH_QUA_4 :
+ case MSH_POLYG_ :
{
std::vector<DI_CuttingPoint> cp;
std::vector<DI_Quad> subQuads;
@@ -605,7 +713,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
T.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
subQuads, subTriangles, boundLines);
}
- else {
+ 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(),
@@ -613,10 +721,21 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
Q.cut(*ls, ipV, ipS, cp, integOrder,integOrder,integOrder,
subQuads, subTriangles, boundLines);
}
+ else if(eType == MSH_POLYG_){
+ for(int i = 0; i < e->getNumChildren(); i++) {
+ MTriangle *t = (MTriangle*) e->getChild(i);
+ DI_Triangle Tri(t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+ t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+ t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z());
+ Tri.cut(*ls, ipV, ipS, cp, integOrder, integOrder, integOrder,
+ subQuads, subTriangles, boundLines);
+ }
+ }
// if cut
if((eType == MSH_TRI_3 && subTriangles.size() > 1) ||
- (eType == MSH_QUA_4 && subTriangles.size() > 2)) {
+ (eType == MSH_QUA_4 && subTriangles.size() > 2) ||
+ (eType == MSH_POLYG_ && subTriangles.size() > e->getNumChildren())) {
std::vector<MTriangle*> poly[2];
for (unsigned int i = 0; i < subTriangles.size(); i++){
@@ -656,7 +775,6 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
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 * -1, 2, physicals);
assignPhysicals(GM, gePhysicals, elementary, 2, physicals);
for (unsigned int i = 0; i < boundLines.size(); i++){
@@ -685,8 +803,11 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
}
}
MLineBorder *lin = new MLineBorder(mv[0], mv[1], p1, p2, ++numEle, ePart);
- border[0][boundLines[i].lsTag()].push_back(lin);
- entityCut[boundLines[i].lsTag()].push_back(ge);
+ 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);
}
}
@@ -696,6 +817,8 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
elements[2][reg].push_back(copy);
else if(eType == MSH_QUA_4)
elements[3][reg].push_back(copy);
+ else if(eType == MSH_POLYG_)
+ elements[8][reg].push_back(copy);
assignPhysicals(GM, gePhysicals, reg, 2, physicals);
}
}
@@ -773,10 +896,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
GModel *cutGM = new GModel(gm->getName() + "_cut");
cutGM->setFileName(cutGM->getName());
- std::map<int, std::vector<MElement*> > border[2];
std::vector<MVertex*> newVertices;
std::vector<GEntity*> gmEntities;
- std::map<int, std::vector<GEntity*> > entityCut;
gm->getEntities(gmEntities);
int numEle = gm->getNumMeshElements();
@@ -784,48 +905,18 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::map<int, int> newtags[4];
for(int i = 0; i < 4; i++)
newtags[i][0] = gm->getMaxElementaryNumber(i);
+ std::map<int, int> borderTags[2];
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, border, physicals, entityCut, newtags);
+ elements, physicals, borderTags, newtags);
cutGM->getMeshPartitions().insert(e->getPartition());
}
}
- // add borders in map elements and change the tag if it's already used
- std::map<int, std::vector<MElement*> >::iterator itbo, itel;
- for(itbo = border[0].begin(); itbo != border[0].end(); itbo++) {
- int reg = itbo->first;
- if(elements[1].count(reg)) {
- itel = elements[1].end(); itel--;
- reg = itel->first + 1;
- }
- for(unsigned int j = 0; j < itbo->second.size(); j++)
- elements[1][reg].push_back(itbo->second[j]);
- std::map<int, std::vector<GEntity*> >::iterator itge = entityCut.find(itbo->first);
- for(unsigned int j = 0; j < itge->second.size(); j++)
- assignPhysicals(gm, itge->second[j]->physicals, reg, 1, physicals);
- }
- for(itbo = border[1].begin(); itbo != border[1].end(); itbo++) {
- int reg = itbo->first;
- if(elements[2].count(reg)) {
- itel = elements[2].end(); itel--;
- reg = itel->first + 1;
- }
- if(elements[3].count(reg)) {
- itel = elements[3].end(); itel--;
- reg = std::max(reg, itel->first + 1);
- }
- for(unsigned int j = 0; j < itbo->second.size(); j++)
- elements[2][reg].push_back(itbo->second[j]);
- std::map<int, std::vector<GEntity*> >::iterator itge = entityCut.find(itbo->first);
- for(unsigned int j = 0; j < itge->second.size(); j++)
- assignPhysicals(gm, itge->second[j]->physicals, reg, 2, physicals);
- }
-
for(unsigned int i = 0; i < newVertices.size(); i++) {
vertexMap[newVertices[i]->getNum()] = newVertices[i];
}
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index e8f302e0b0cc00eba0c53fd031b206d43997a118..bd262edc5653d5433a85513dc85b0565f1cf3b09 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -263,6 +263,28 @@ class MTriangleBorder : public MTriangle {
}
~MTriangleBorder() {}
MPolyhedron* getDomain(int i) const { return _domains[i]; }
+ virtual MElement *getParent() const { return _domains[0]->getParent(); }
+ virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
+ {
+ return getParent()->getFunctionSpace(order);
+ }
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+ {
+ getParent()->getShapeFunctions(u, v, w, s, o);
+ }
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+ {
+ getParent()->getGradShapeFunctions(u, v, w, s, o);
+ }
+ virtual void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ getParent()->xyz2uvw(xyz,uvw);
+ }
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual int getNumIntegrationPointsToAllocate (int pOrder)
+ {
+ return getNGQTPts(pOrder);
+ }
};
class MLineBorder : public MLine {
@@ -277,6 +299,28 @@ class MLineBorder : public MLine {
}
~MLineBorder() {}
MPolygon* getDomain(int i) const { return _domains[i]; }
+ virtual MElement *getParent() const { return _domains[0]->getParent(); }
+ virtual const gmshFunctionSpace* getFunctionSpace(int order=-1) const
+ {
+ return getParent()->getFunctionSpace(order);
+ }
+ virtual void getShapeFunctions(double u, double v, double w, double s[], int o)
+ {
+ getParent()->getShapeFunctions(u, v, w, s, o);
+ }
+ virtual void getGradShapeFunctions(double u, double v, double w, double s[][3], int o)
+ {
+ getParent()->getGradShapeFunctions(u, v, w, s, o);
+ }
+ virtual void xyz2uvw(double xyz[3], double uvw[3])
+ {
+ getParent()->xyz2uvw(xyz,uvw);
+ }
+ virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const;
+ virtual int getNumIntegrationPointsToAllocate (int pOrder)
+ {
+ return pOrder / 2 + 1;
+ }
};
GModel *buildCutMesh(GModel *gm, gLevelset *ls,
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 8f4689f5619befa55121ef5787d2812631c5a8aa..9b7c4b0dc5620520a0b4ac8275907586652c08fe 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -37,6 +37,11 @@ class MLine : public MElement {
virtual int getDim(){ return 1; }
virtual int getNumVertices() const { return 2; }
virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual void getCoordinates (int num, double c[3]) const {
+ c[0] = _v[num]->x();
+ c[1] = _v[num]->y();
+ c[2] = _v[num]->z();
+ }
virtual int getNumEdges(){ return 1; }
virtual MEdge getEdge(int num){ return MEdge(_v[0], _v[1]); }
virtual int getNumEdgesRep(){ return 1; }
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 13516f421cce53a15d59dc141f0d0c04541f3fa4..d5fbe71c556407d75cadd19d800f43e2d64cee3c 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -54,6 +54,11 @@ class MTriangle : public MElement {
virtual double distoShapeMeasure();
virtual int getNumVertices() const { return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
+ virtual void getCoordinates (int num, double c[3]) const {
+ c[0] = _v[num]->x();
+ c[1] = _v[num]->y();
+ c[2] = _v[num]->z();
+ }
virtual MVertex *getVertexMED(int num)
{
static const int map[3] = {0, 2, 1};
diff --git a/contrib/DiscreteIntegration/DILevelset.h b/contrib/DiscreteIntegration/DILevelset.h
index c4bfa4072041111be36ae5468365f789a8f38d20..2647cbb72745db146c19bcc61c3ce836c13372a8 100644
--- a/contrib/DiscreteIntegration/DILevelset.h
+++ b/contrib/DiscreteIntegration/DILevelset.h
@@ -3,7 +3,6 @@
#include <cmath>
#include <vector>
-#include <cassert>
#include <stdio.h>
@@ -80,7 +79,7 @@ public:
gLevelsetPrimitive(int &tag) {
if (tag < 1) {
printf("Tag of the levelset (%d) must be greater than 0.\n", tag);
- throw;
+ tag = std::fabs(tag);
}
tag_ = tag++;
}
@@ -88,7 +87,7 @@ public:
std::vector<const gLevelset *> getChildren() const { std::vector<const gLevelset *> p; return p; }
double choose (double d1, double d2) const {
printf("Cannot use function \"choose\" with a primitive!\n");
- throw;
+ return d1;
}
virtual int type() const = 0;
bool isPrimitive() const {return true;}
diff --git a/contrib/DiscreteIntegration/Integration3D.cpp b/contrib/DiscreteIntegration/Integration3D.cpp
index 2b18696eb4911cbf9ad03587912e79c159ce1f43..ad888429c3c54586327614ae862bc19b61d69dc7 100644
--- a/contrib/DiscreteIntegration/Integration3D.cpp
+++ b/contrib/DiscreteIntegration/Integration3D.cpp
@@ -184,7 +184,7 @@ int commonV (int &s11, int &s12, int &s21, int &s22) {
if(s11 == s21 || s11 == s22) return s11;
if(s12 == s21 || s12 == s22) return s12;
printf("no common summit, %d,%d,%d,%d\n", s11, s12, s21, s22);
- throw;
+ return 0;
}
double adjustLs (double ls) {
@@ -218,7 +218,6 @@ int minimum(double *x, double *y, double *z, const int num) {
for(int i = 1; i < county; i++) if(z[INDy[i]] < zm) zm = z[INDy[i]];
std::vector<int> INDz(county); int countz = 0;
for(int i = 0; i < county; i++) if(z[INDy[i]] == zm) INDz[countz++] = INDy[i];
- assert (countz == 1);
return INDz[0];
}
@@ -410,38 +409,37 @@ int bestQuality (const DI_Point &p0, const DI_Point &p1, const DI_Point &p2,
}
}
- assert(cut > -1);
if(cut == 0) {
t1 = DI_Tetra(p0, p1, p2, p5);
t2 = DI_Tetra(p0, p1, p5, p3);
t3 = DI_Tetra(p1, p5, p3, p4);
return 1;
}
- if(cut == 1) {
+ else if(cut == 1) {
t1 = DI_Tetra(p0, p1, p2, p5);
t2 = DI_Tetra(p0, p1, p5, p4);
t3 = DI_Tetra(p0, p4, p5, p3);
return 2;
}
- if(cut == 2) {
+ else if(cut == 2) {
t1 = DI_Tetra(p0, p1, p2, p4);
t2 = DI_Tetra(p0, p4, p2, p5);
t3 = DI_Tetra(p0, p4, p5, p3);
return 3;
}
- if(cut == 3) {
+ else if(cut == 3) {
t1 = DI_Tetra(p0, p1, p2, p4);
t2 = DI_Tetra(p0, p4, p2, p3);
t3 = DI_Tetra(p2, p3, p4, p5);
return 4;
}
- if(cut == 4) {
+ else if(cut == 4) {
t1 = DI_Tetra(p0, p1, p2, p3);
t2 = DI_Tetra(p1, p2, p3, p4);
t3 = DI_Tetra(p2, p3, p4, p5);
return 5;
}
- else { //cut == 5
+ else if(cut == 5) {
t1 = DI_Tetra(p0, p1, p2, p3);
t2 = DI_Tetra(p1, p2, p3, p5);
t3 = DI_Tetra(p1, p5, p3, p4);
@@ -480,9 +478,11 @@ DI_Point quadMidNode(const DI_Point &p1, const DI_Point &p2, const DI_Point &pf,
double v1f[3]; vec(p1, pf, v1f);
double nf[3]; cross(v12, v1f, nf);
double bisector[3]; cross(nf, v12, bisector);
- double normB = norm(bisector); assert (normB != 0);
- for (int i = 0; i < 3; i++)
- bisector[i] = bisector[i] / normB;
+ double normB = norm(bisector);
+ if(normB) {
+ for (int i = 0; i < 3; i++)
+ bisector[i] = bisector[i] / normB;
+ }
// raise the length of bisector if needed .........
DI_Point pt(midEN.x() + bisector[0], midEN.y() + bisector[1], midEN.z() + bisector[2]);
pt.addLs(e);
@@ -548,13 +548,12 @@ void DI_IntegrationPoint::computeLs (const DI_Element *e, const std::vector<cons
Ls.push_back(ls);
}
}
- assert (Ls.size() == 1);
setLs(Ls.back());
}
// DI_CuttingPoint methods -----------------------------------------------------------------------------------------------------
void DI_CuttingPoint::chooseLs (const gLevelset *Lsi) {
- assert(Ls.size() > 1);
+ if(Ls.size() < 2) return;
double ls1 = Ls[Ls.size() - 2], ls2 = Ls[Ls.size() - 1];
double ls = Lsi->choose(ls1, ls2);
Ls.pop_back(); Ls.pop_back();
@@ -588,7 +587,10 @@ DI_Element::DI_Element(const DI_Element &cp) : lsTag_(cp.lsTag()), polOrder_(cp.
if(mid_) delete mid_;
}*/
DI_Element & DI_Element::operator= (const DI_Element &rhs){
- assert (type() == rhs.type());
+ if(type() != rhs.type()) {
+ printf("Error : try to assign element of different type!\n");
+ return *this;
+ }
if(this != &rhs) {
for(int i = 0; i < nbVert(); i++) {
delete pts_[i];
@@ -659,14 +661,16 @@ void DI_Element::addLs (const double *ls) {
mid_[i]->addLs(ls[nbVert()+i]);
}
void DI_Element::addLs (const DI_Element *e) {
- assert(e->sizeLs() > 0);
+ if(e->sizeLs() < 1) return;
for(int i = 0; i < nbVert(); i++)
pts_[i]->addLs(e);
for(int i = 0; i < nbMid(); i++)
mid_[i]->addLs(e);
}
void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls) {
- assert (type() == e->type());
+ if(type() != e->type()) {
+ printf("Error : addLs with element of different type\n");
+ }
for(int j = 0; j < nbVert(); ++j) {
double ls = Ls(e->x(j), e->y(j), e->z(j));
pts_[j]->addLs(adjustLs(ls));
@@ -683,7 +687,7 @@ void DI_Element::addLs (const DI_Element *e, const gLevelset &Ls) {
}
void DI_Element::chooseLs (const gLevelset *Lsi) {
if(sizeLs() < 2)
- printf("chooseLs with element size < 2 : typeEl=%d\n", type());
+ printf("chooseLs with element ls size < 2 : typeEl=%d\n", type());
for(int i = 0; i < nbVert(); i++)
pts_[i]->chooseLs(Lsi);
for(int i = 0; i < nbMid(); i++)
@@ -696,7 +700,7 @@ void DI_Element::clearLs() {
(mid_[i]->Ls).clear();
}
bool DI_Element::addQuadEdge (int edge, DI_Point *xm, const DI_Element *e, const std::vector<const gLevelset *> RPNi) {
- /*if(edge >= nbEdg()) {printf("wrong number (%d) for quadratic edge for a ", edge); print(); throw;}
+ /*if(edge >= nbEdg()) {printf("wrong number (%d) for quadratic edge for a ", edge); print(); return false;}
int s1, s2; vert(edge, s1, s2);
bool quad0 = isQuad();
if(!quad0) quad(e, RPNi);
@@ -753,7 +757,7 @@ bool DI_Element::contain (const DI_Point &p) const {
}
return true;
default :
- throw;
+ return false;
}
}
bool DI_Element::contain (const DI_Element *e) const {
@@ -912,10 +916,10 @@ void DI_Element::computeLsTagDom(const DI_Element *e, const std::vector<const gL
if(mid.isInsideDomain())
{setLsTag(1); return;}
}
- printf("Unable to determine the sign of the element : \n");
+ printf("Error : Unable to determine the sign of the element : \n");
printf("Parent element : "); e->printls();
printf("Element : "); printls();
- throw;
+ return;
}
// set the lsTag to -1 if the element is not on the border of the domain
void DI_Element::computeLsTagBound(std::vector<DI_Hexa> &hexas, std::vector<DI_Tetra> &tetras){
@@ -1611,7 +1615,8 @@ void DI_Triangle::selfSplit (const DI_Element *e, const std::vector<const gLevel
break;
}
default:
- throw;
+ printf("Error : %d edge(s) cut in the triangle (ls : %g %g %g)\n",
+ COUNT, pt(0).ls(), pt(1).ls(), pt(2).ls());
}
}
@@ -1963,8 +1968,8 @@ void DI_Tetra::selfSplit (const DI_Element *e, const std::vector<const gLevelset
break;
}
default:
- printf("Error : %d edge(s) cut (ls : %g %g %g %g)\n", COUNT, pt(0).ls(), pt(1).ls(), pt(2).ls(), pt(3).ls());
- throw;
+ printf("Error : %d edge(s) cut in the tetrahedron (ls : %g %g %g %g)\n",
+ COUNT, pt(0).ls(), pt(1).ls(), pt(2).ls(), pt(3).ls());
}
}
@@ -2091,7 +2096,6 @@ void DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
}
for(int l = 0; l < (int)ll_subLines.size(); l++) {
- assert(ll_subLines[l].sizeLs() == 1);
ll_subLines[l].computeLsTagDom(&ll, RPN);
DI_Line ll_subLn = ll_subLines[l];
mappingEl(&ll_subLn);
@@ -2099,7 +2103,6 @@ void DI_Line::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
lines.push_back(ll_subLn);
}
for(int p = 0; p < (int)ll_cp.size(); p++) {
- assert (ll_cp[p].sizeLs() == 1);
if(ll_cp[p].ls() != 0) continue;
mappingCP(ll_cp[p]);
bool isIn = false;
@@ -2265,7 +2268,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++) {
- assert(tt_subQuads[q].sizeLs() == 1);
tt_subQuads[q].computeLsTagDom(&tt, RPN);
DI_Quad tt_subQ = tt_subQuads[q];
mappingEl(&tt_subQ);
@@ -2273,7 +2275,6 @@ void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
subQuads.push_back(tt_subQ);
}
for(int t = 0; t < (int)tt_subTriangles.size(); t++) {
- assert(tt_subTriangles[t].sizeLs() == 1);
tt_subTriangles[t].computeLsTagDom(&tt, RPN);
DI_Triangle tt_subTr = tt_subTriangles[t];
mappingEl(&tt_subTr);
@@ -2281,7 +2282,6 @@ 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++) {
- assert(tt_surfLines[l].sizeLs() == 1);
tt_surfLines[l].computeLsTagBound(tt_subQuads, tt_subTriangles); //tt_surfLines[l].printls();
if(tt_surfLines[l].lsTag() == -1) continue;
DI_Line tt_surfLn = tt_surfLines[l];
@@ -2290,7 +2290,6 @@ void DI_Triangle::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip
surfLines.push_back(tt_surfLn);
}
for(int p = 0; p < (int)tt_cp.size(); p++) {
- assert (tt_cp[p].sizeLs() == 1);
if(tt_cp[p].ls() != 0) continue;
mappingCP(tt_cp[p]);
bool isIn = false;
@@ -2468,7 +2467,6 @@ void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
}*/
for(int q = 0; q < (int)qq_subQuads.size(); q++) {
- assert(qq_subQuads[q].sizeLs() == 1);
qq_subQuads[q].computeLsTagDom(&qq, RPN);
DI_Quad qq_subQ = qq_subQuads[q];
mappingEl(&qq_subQ);
@@ -2476,7 +2474,6 @@ void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
subQuads.push_back(qq_subQ);
}
for(int t = 0; t < (int)qq_subTriangles.size(); t++) {
- assert(qq_subTriangles[t].sizeLs() == 1);
qq_subTriangles[t].computeLsTagDom(&qq, RPN);
DI_Triangle qq_subTr = qq_subTriangles[t];
mappingEl(&qq_subTr); //qq_subTr.printls();
@@ -2484,7 +2481,6 @@ void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
subTriangles.push_back(qq_subTr);
}
for(int l = 0; l < (int)qq_surfLines.size(); l++) {
- assert(qq_surfLines[l].sizeLs() == 1);
qq_surfLines[l].computeLsTagBound(qq_subQuads, qq_subTriangles);
if(qq_surfLines[l].lsTag() == -1) continue; //FIXME
DI_Line qq_surfLn = qq_surfLines[l];
@@ -2493,7 +2489,6 @@ void DI_Quad::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
surfLines.push_back(qq_surfLn);
}
for(int p = 0; p < (int)qq_cp.size(); p++) {
- assert (qq_cp[p].sizeLs() == 1);
if(qq_cp[p].ls() != 0) continue;
mappingCP(qq_cp[p]);
bool isIn = false;
@@ -2660,7 +2655,6 @@ void DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
QError[i].print(this);
for(int t = 0; t < (int)tt_subTetras.size(); t++) {
- assert(tt_subTetras[t].sizeLs() == 1);
tt_subTetras[t].computeLsTagDom(&tt, RPN);
DI_Tetra tt_subT = tt_subTetras[t];
mappingEl(&tt_subT);
@@ -2668,7 +2662,6 @@ void DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
subTetras.push_back(tt_subT);
}
for(int q = 0; q < (int)tt_surfQuads.size(); q++) {
- assert(tt_surfQuads[q].sizeLs() == 1);
tt_surfQuads[q].computeLsTagBound(tt_subHexas, tt_subTetras);
if(tt_surfQuads[q].lsTag() == -1) continue;
DI_Quad tt_surfQ = tt_surfQuads[q];
@@ -2677,7 +2670,6 @@ void DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
surfQuads.push_back(tt_surfQ);
}
for(int t = 0; t < (int)tt_surfTriangles.size(); t++) {
- assert(tt_surfTriangles[t].sizeLs() == 1);
tt_surfTriangles[t].computeLsTagBound(tt_subHexas, tt_subTetras);
if(tt_surfTriangles[t].lsTag() == -1) continue;
DI_Triangle tt_surfTr = tt_surfTriangles[t];
@@ -2686,7 +2678,6 @@ void DI_Tetra::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
surfTriangles.push_back(tt_surfTriangles[t]);
}
for(int p = 0; p < (int)tt_cp.size(); p++) {
- assert (tt_cp[p].sizeLs() == 1);
if(tt_cp[p].ls() != 0) continue;
mappingCP(tt_cp[p]);
bool isIn = false;
@@ -2901,7 +2892,6 @@ void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
QError[i].print(this);
for(int h = 0; h < (int)hh_subHexas.size(); h++) {
- assert (hh_subHexas[h].sizeLs() == 1);
hh_subHexas[h].computeLsTagDom(&hh, RPN);
DI_Hexa hh_subH = hh_subHexas[h];
mappingEl(&hh_subH);
@@ -2909,7 +2899,6 @@ void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
subHexas.push_back(hh_subH);
}
for(int t = 0; t < (int)hh_subTetras.size(); t++) {
- assert (hh_subTetras[t].sizeLs() == 1);
hh_subTetras[t].computeLsTagDom(&hh, RPN);
DI_Tetra hh_subT = hh_subTetras[t];
mappingEl(&hh_subT);
@@ -2917,7 +2906,6 @@ void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
subTetras.push_back(hh_subT);
}
for(int q = 0; q < (int)hh_surfQuads.size(); q++) {
- assert(hh_surfQuads[q].sizeLs() == 1);
hh_surfQuads[q].computeLsTagBound(hh_subHexas, hh_subTetras);
if(hh_surfQuads[q].lsTag() == -1) continue;
DI_Quad hh_surfQ = hh_surfQuads[q];
@@ -2926,7 +2914,6 @@ void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
surfQuads.push_back(hh_surfQ);
}
for(int t = 0; t < (int)hh_surfTriangles.size(); t++) {
- assert (hh_surfTriangles[t].sizeLs() == 1);
hh_surfTriangles[t].computeLsTagBound(hh_subHexas, hh_subTetras);
if(hh_surfTriangles[t].lsTag() == -1) continue;
DI_Triangle hh_surfTr = hh_surfTriangles[t];
@@ -2935,7 +2922,6 @@ void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
surfTriangles.push_back(hh_surfTr);
}
for(int l = 0; l < (int)hh_frontLines.size(); l++) {
- assert(hh_frontLines[l].sizeLs() == 1);
hh_frontLines[l].computeLsTagBound(hh_surfQuads, hh_surfTriangles);
if(hh_frontLines[l].lsTag() == -1) continue;
DI_Line hh_frontLn = hh_frontLines[l];
@@ -2944,7 +2930,6 @@ void DI_Hexa::cut (const gLevelset &Ls, std::vector<DI_IntegrationPoint> &ip,
frontLines.push_back(hh_frontLn);
}
for(int p = 0; p < (int)hh_cp.size(); p++) {
- assert (hh_cp[p].sizeLs() == 1);
if(hh_cp[p].ls() != 0) continue; // returns only the cutting points with ls==0
mappingCP(hh_cp[p]);
bool isIn = false;
@@ -2986,7 +2971,8 @@ void DI_Hexa::cut (const DI_Element *e, const std::vector<const gLevelset *> RPN
else{
if(on == 4){ // the level set is zero on a face of the hex
// assert the nodes are in the same plane
- if(!isPlanar(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) { printf("THE FOUR NODES WITH ZERO LS ARE NOT PLANAR"); throw;}
+ if(!isPlanar(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) {
+ printf("Error : The 4 nodes with zero levelset are not planar!\n"); return;}
// order the 4 nodes
if(!ordered4Nodes(pt(ze[0]), pt(ze[1]), pt(ze[2]), pt(ze[3]))) swap(ze[2], ze[3]);
// add the quad twice if the face belongs to 2 hexas => remove it later!
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 1e55ac6ba0cb8c48082a6d48c0eb0209ec3966af..862e040ab497a96333a1cc1b74638cb9a3b6296d 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -3,7 +3,6 @@
#include <list>
#include <vector>
-#include <cassert>
#include <cmath>
#include "DILevelset.h"
@@ -164,7 +163,7 @@ static inline double TetraVol(double x1, double y1, double z1, double x2, double
double vol = ((x2 - x1) * ((y3 - y1) * (z4 - z1) - (y4 - y1) * (z3 - z1))
- (x3 - x1) * ((y2 - y1) * (z4 - z1) - (y4 - y1) * (z2 - z1))
+ (x4 - x1) * ((y2 - y1) * (z3 - z1) - (y3 - y1) * (z2 - z1))) / 6.;
- if(vol < 0) {printf("TET HAS NEGATIVE VOLUME = %g\n", vol); throw;}
+ if(vol < 0) {printf("TET HAS NEGATIVE VOLUME = %g\n", vol);}
return vol;
}
static inline double TetraVol(const DI_Point p1, const DI_Point p2,
@@ -849,7 +848,8 @@ public:
if(i == 1) return p21_;
if(i == 2) return p12_;
if(i == 3) return p22_;
- throw;
+ printf("DI_QualError::pt only accept indices from 0 to 3!\n");
+ DI_Point p; return p;
}
void print(const DI_Element *e) const{
DI_Point pt1 = p11_, pt2 = p12_, pt3 = p21_, pt4 = p22_;