diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 8d61a2ee2b34c5217d7c9a276010abba0d05ed57..7b663f55595dd710eab810089157e7d0fb2e65cc 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -3,12 +3,11 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to <gmsh@geuz.org>.
+#include "GmshConfig.h"
#include "GModel.h"
#include "MElement.h"
#include "MElementCut.h"
-//#define HAVE_DINTEGRATION
-
#if defined(HAVE_DINTEGRATION)
#include "DILevelset.h"
#include "Integration3D.h"
@@ -72,6 +71,7 @@ void MPolyhedron::_init()
}
}
}
+
bool MPolyhedron::isInside(double u, double v, double w)
{
double ksi[3] = {u, v, w};
@@ -83,6 +83,7 @@ bool MPolyhedron::isInside(double u, double v, double w)
}
return false;
}
+
void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
{
*npts = 0;
@@ -93,7 +94,8 @@ void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
_parts[i]->getIntegrationPoints(pOrder, &nptsi, &ptsi);
double uvw[4][3];
for(int j = 0; j < 4; j++) {
- double xyz[3] = {_parts[i]->getVertex(j)->x(), _parts[i]->getVertex(j)->y(),
+ double xyz[3] = {_parts[i]->getVertex(j)->x(),
+ _parts[i]->getVertex(j)->y(),
_parts[i]->getVertex(j)->z()};
_orig->xyz2uvw(xyz, uvw[j]);
}
@@ -118,8 +120,9 @@ void MPolyhedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
*npts += nptsi;
}
}
+
void MPolyhedron::writeMSH(FILE *fp, double version, bool binary, int num,
- int elementary, int physical)
+ int elementary, int physical)
{
int type = getTypeForMSH();
@@ -287,6 +290,7 @@ bool MPolygon::isInside(double u, double v, double w)
}
return false;
}
+
void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
{
*npts = 0;
@@ -320,6 +324,7 @@ void MPolygon::getIntegrationPoints(int pOrder, int *npts, IntPt **pts) const
*npts += nptsi;
}
}
+
void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
int elementary, int physical)
{
@@ -370,18 +375,23 @@ void MPolygon::writeMSH(FILE *fp, double version, bool binary, int num,
#if defined(HAVE_DINTEGRATION)
-int getElementVertexNum (DI_Point p, MElement *e)
+static bool equalV(MVertex *v, DI_Point p)
+{
+ return (fabs(v->x() - p.x()) < 1.e-15 &&
+ fabs(v->y() - p.y()) < 1.e-15 &&
+ fabs(v->z() - p.z()) < 1.e-15);
+}
+
+static int getElementVertexNum(DI_Point p, MElement *e)
{
for(int i = 0; i < e->getNumVertices(); i++)
- if(fabs(p.x() - e->getVertex(i)->x()) < 1.e-12 &&
- fabs(p.y() - e->getVertex(i)->y()) < 1.e-12 &&
- fabs(p.z() - e->getVertex(i)->z()) < 1.e-12)
+ if(equalV(e->getVertex(i), p))
return e->getVertex(i)->getNum();
return -1;
}
-void assignPhysicals(GModel *GM, std::vector<int> gePhysicals, int reg, int dim,
- std::map<int, std::map<int, std::string> > physicals[4])
+static void assignPhysicals(GModel *GM, std::vector<int> gePhysicals, int reg, int dim,
+ std::map<int, std::map<int, std::string> > physicals[4])
{
for(unsigned int i = 0; i < gePhysicals.size(); i++){
int phys = gePhysicals[i];
@@ -390,12 +400,6 @@ void assignPhysicals(GModel *GM, std::vector<int> gePhysicals, int reg, int dim,
}
}
-bool equalV(MVertex *v, DI_Point p)
-{
- return (fabs(v->x() - p.x()) < 1.e-15 && fabs(v->y() - p.y()) < 1.e-15 &&
- fabs(v->z() - p.z()) < 1.e-15);
-}
-
static int getElementaryTag(double ls, int elementary, std::map<int, int> &newtags)
{
if(ls < 0){
@@ -775,7 +779,7 @@ static void elementCutMesh(MElement *e, gLevelset *ls, GEntity *ge, GModel *GM,
break;
default :
Msg::Error("This type of element cannot be cut.");
- throw;
+ return;
}
}
@@ -788,6 +792,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
{
#if defined(HAVE_DINTEGRATION)
GModel *cutGM = new GModel(gm->getName() + "_cut");
+ cutGM->setFileName(cutGM->getName());
std::map<int, std::vector<MElement*> > border[2];
std::vector<MVertex*> newVertices;
@@ -842,11 +847,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
assignPhysicals(gm, itge->second[j]->physicals, reg, 2, physicals);
}
- // number the new vertices and add in vertexMap
- //std::map<int, MVertex* >::iterator itend = vertexMap.end(); itend--;
- //int num = itend->first;
for(unsigned int i = 0; i < newVertices.size(); i++) {
- //newVertices[i]->setNum(++num);
vertexMap[newVertices[i]->getNum()] = newVertices[i];
}
printf("numbering vertices finished : %d vertices \n", (int)vertexMap.size());
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index fe7d524010432171e5ef3f72c3b358256695cd8d..65b594f48fa6da981530c10cc3a3b3bc5618a937 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -6,6 +6,7 @@
#ifndef _MELEMENTCUT_H_
#define _MELEMENTCUT_H_
+#include "GmshMessage.h"
#include "MElement.h"
#include "MTetrahedron.h"
#include "MTriangle.h"
@@ -26,7 +27,16 @@ class MPolyhedron : public MElement {
void _init();
public:
MPolyhedron(std::vector<MVertex*> v, int num=0, int part=0)
- : MElement(num, part), _owner(false), _orig(0) {}
+ : MElement(num, part), _owner(false), _orig(0)
+ {
+ if(v.size() % 4){
+ Msg::Error("Got %d vertices for polyhedron", (int)v.size());
+ return;
+ }
+ for(unsigned int i = 0; i < v.size(); i += 4)
+ _parts.push_back(new MTetrahedron(v[i], v[i + 1], v[i + 2], v[i + 3]));
+ _init();
+ }
MPolyhedron(std::vector<MTetrahedron*> vT, int num=0, int part=0)
: MElement(num, part), _owner(false), _orig(0)
{
diff --git a/benchmarks/levelset/carreTri.geo b/benchmarks/levelset/carreTri.geo
index d05a7c6bd3a438282921b4787be66fb29c58d5ea..69228ab113a9d3c5741771c5eac3230fcbfc6b49 100644
--- a/benchmarks/levelset/carreTri.geo
+++ b/benchmarks/levelset/carreTri.geo
@@ -19,10 +19,10 @@ Line Loop(5) = {1,2,3,4} ;
Plane Surface(6) = {5} ;
//Recombine Surface{6};
-//Physical Line(100) = {1};
-//Physical Line(200) = {2};
-//Physical Line(300) = {3};
-//Physical Line(400) = {4};
+Physical Line(100) = {1};
+Physical Line(200) = {2};
+Physical Line(300) = {3};
+Physical Line(400) = {4};
//Physical Surface(1000) = {6};
Mesh 2;