From 5dd35e2aef46261d6e93e340c93ceed21eee93d7 Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Wed, 10 Feb 2010 10:50:20 +0000
Subject: [PATCH] SplitMesh
---
Geo/GModel.cpp | 4 +-
Geo/GModel.h | 3 +-
Geo/GModelIO_Mesh.cpp | 6 +-
Geo/MElementCut.cpp | 197 ++++++++++++++++----
Geo/MElementCut.h | 3 +-
Geo/MVertex.h | 1 +
Numeric/fullMatrix.h | 4 +-
Parser/Gmsh.y | 6 +
contrib/DiscreteIntegration/Integration3D.h | 3 -
9 files changed, 182 insertions(+), 45 deletions(-)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index ac28226d9d..0956f1ee77 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1302,7 +1302,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
}
}
-GModel *GModel::buildCutGModel(gLevelset *ls)
+GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem)
{
std::map<int, std::vector<MElement*> > elements[10];
std::map<int, std::map<int, std::string> > physicals[4];
@@ -1311,7 +1311,7 @@ GModel *GModel::buildCutGModel(gLevelset *ls)
Msg::Info("Cutting mesh...");
double t1 = Cpu();
- GModel *cutGM = buildCutMesh(this, ls, elements, vertexMap, physicals);
+ GModel *cutGM = buildCutMesh(this, ls, elements, vertexMap, physicals, cutElem);
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
cutGM->_storeElementsInEntities(elements[i]);
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 9d7d17a266..af292e094b 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -340,7 +340,8 @@ class GModel
int mesh(int dimension);
// build a new GModel by cutting the elements crossed by the levelset ls
- GModel *buildCutGModel(gLevelset *ls);
+ // if cutElem is set to false, split the model without cutting the elements
+ GModel *buildCutGModel(gLevelset *ls, bool cutElem = true);
// create a GModel by importing a mesh (vertexMap has a dim equal to
// the number of vertices and all the other vectors have a dim equal
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 811e843869..e569fd11d1 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -318,7 +318,7 @@ int GModel::readMSH(const std::string &name)
if(numVertices > 100000)
Msg::ProgressMeter(i + 1, numVertices, "Reading nodes");
}
- // if the vertex numbering is dense, tranfer the map into a
+ // If the vertex numbering is dense, transfer the map into a
// vector to speed up element creation
if((int)vertexMap.size() == numVertices &&
((minVertex == 1 && maxVertex == numVertices) ||
@@ -2823,7 +2823,7 @@ int GModel::readDIFF(const std::string &name)
vertexMap[num] = new MVertex(xyz[0], xyz[1], xyz[2], 0, num);
if(numVertices > 100000)
Msg::ProgressMeter(i + 1, numVertices, "Reading nodes");
- // If the vertex numbering is dense, tranfer the map into a
+ // If the vertex numbering is dense, transfer the map into a
// vector to speed up element creation
if((int)vertexMap.size() == numVertices &&
((minVertex == 1 && maxVertex == numVertices) ||
@@ -3175,7 +3175,7 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
if(minVertex == std::numeric_limits<int>::max())
Msg::Error("Could not determine the min index of vertices");
- // If the vertex numbering is dense, tranfer the map into a
+ // If the vertex numbering is dense, transfer the map into a
// vector to speed up element creation
if((minVertex == 1 && maxVertex == numVertices) ||
(minVertex == 0 && maxVertex == numVertices - 1)){
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index f82d2ba7f4..357cce3614 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -384,8 +384,6 @@ void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
//---------------------------------------- CutMesh ----------------------------
-#if defined(HAVE_DINTEGRATION)
-
static bool equalV(MVertex *v, const DI_Point *p)
{
return (fabs(v->x() - p->x()) < 1.e-15 &&
@@ -424,7 +422,8 @@ static int getElementaryTag(int lsTag, int elementary, std::map<int, int> &newEl
}
return elementary;
}
-static void getPhysicalTag(int lsTag, const std::vector<int> &physicals, std::vector<int> &phys2, std::map<int, int> &newPhysTags)
+static void getPhysicalTag(int lsTag, const std::vector<int> &physicals,
+ std::vector<int> &phys2, std::map<int, int> &newPhysTags)
{
phys2.clear();
for(unsigned int i = 0; i < physicals.size(); i++){
@@ -452,11 +451,127 @@ static int getBorderTag(int lsTag, int count, int &maxTag, std::map<int, int> &b
return lsTag;
}
-typedef std::set<MVertex*,MVertexLessThanLexicographic> newVerticesContainer ;
+static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
+ GEntity *ge, GModel *GM, int &numEle,
+ std::map<int, MVertex*> &vertexMap,
+ std::map<int, std::vector<MElement*> > elements[10],
+ std::map<int, std::map<int, std::string> > physicals[4],
+ std::map<int, int> newElemTags[4],
+ std::map<int, int> newPhysTags[4])
+{
+ int elementary = ge->tag();
+ int eType = e->getTypeForMSH();
+ int ePart = e->getPartition();
+ std::vector<int> gePhysicals = ge->physicals;
+
+ std::vector<MVertex*> vmv;
+ 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;
+ MElement *copy = factory.create(eType, vmv, ++numEle, ePart);
+
+ double lsMean = 0.;
+ for(int k = 0; k < e->getNumVertices(); k++)
+ lsMean += verticesLs(e->getVertex(k)->getIndex(), 0);
+ int lsTag = (lsMean < 0) ? 1 : -1;
+
+ switch (eType) {
+ case MSH_TET_4 :
+ case MSH_HEX_8 :
+ case MSH_PYR_5 :
+ case MSH_PRI_6 :
+ case MSH_POLYH_ :
+ {
+ int reg = getElementaryTag(lsTag, elementary, newElemTags[3]);
+ std::vector<int> phys;
+ getPhysicalTag(lsTag, gePhysicals, phys, newPhysTags[3]);
+ if(eType == MSH_TET_4)
+ elements[4][reg].push_back(copy);
+ else if(eType == MSH_HEX_8)
+ elements[5][reg].push_back(copy);
+ else if(eType == MSH_PRI_6)
+ 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, phys, reg, 3, physicals);
+ }
+ break;
+ case MSH_TRI_3 :
+ case MSH_QUA_4 :
+ case MSH_POLYG_ :
+ {
+ int reg = getElementaryTag(lsTag, elementary, newElemTags[2]);
+ std::vector<int> phys;
+ getPhysicalTag(lsTag, gePhysicals, phys, newPhysTags[2]);
+ if(eType == MSH_TRI_3)
+ 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, phys, reg, 2, physicals);
+ }
+ break;
+ case MSH_LIN_2 :
+ {
+ int reg = getElementaryTag(lsTag, elementary, newElemTags[1]);
+ std::vector<int> phys;
+ getPhysicalTag(lsTag, gePhysicals, phys, newPhysTags[1]);
+ elements[1][reg].push_back(copy);
+ assignPhysicals(GM, phys, reg, 1, physicals);
+ }
+ break;
+ case MSH_PNT :
+ {
+ int reg = getElementaryTag(lsTag, elementary, newElemTags[0]);
+ std::vector<int> phys;
+ getPhysicalTag(lsTag, gePhysicals, phys, newPhysTags[0]);
+ elements[0][reg].push_back(copy);
+ assignPhysicals(GM, phys, reg, 0, physicals);
+ }
+ break;
+ default :
+ Msg::Error("This type of element cannot be cut.");
+ return;
+ }
+}
+
+#if defined(HAVE_DINTEGRATION)
+
+typedef std::set<MVertex*, MVertexLessThanLexicographic> newVerticesContainer ;
static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
fullMatrix<double> &verticesLs,
- GEntity *ge, GModel *GM, int &numEle, std::map<int, MVertex*> &vertexMap,
+ GEntity *ge, GModel *GM, int &numEle,
+ std::map<int, MVertex*> &vertexMap,
newVerticesContainer &newVertices,
std::map<int, std::vector<MElement*> > elements[10],
std::map<int, std::map<int, std::string> > physicals[4],
@@ -620,10 +735,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
for(int j = 0; j < 4; j++){
int numV = getElementVertexNum(tetras[i]->pt(j), e);
if (numV == -1){
- MVertex *newv = new MVertex(tetras[i]->x(j), tetras[i]->y(j), tetras[i]->z(j), 0);
+ MVertex *newv = new MVertex(tetras[i]->x(j), tetras[i]->y(j), tetras[i]->z(j));
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
- if (!it.second) delete newv;
+ if (!it.second) newv->deleteLast();
}
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
@@ -677,10 +792,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
for(int j = 0; j < 3; j++){
int numV = getElementVertexNum(triangles[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j), triangles[i]->z(j), 0);
+ MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j), triangles[i]->z(j));
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
- if (!it.second) delete newv;
+ if (!it.second) newv->deleteLast();
}
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
@@ -696,8 +811,8 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
if(p1 || p2) tri = new MTriangleBorder(mv[0], mv[1], mv[2], p1, p2, ++numEle, ePart);
else tri = new MTriangle(mv[0], mv[1], mv[2], ++numEle, ePart);
int lsT = triangles[i]->lsTag();
- int c = elements[2].count(lsT) + elements[3].count(lsT);
- // suppose that the surfaces have been cut before the volumes!
+ int c = elements[2].count(lsT) + elements[3].count(lsT) + elements[8].count(lsT);
+ // the surfaces are cut before the volumes!
int reg = getBorderTag(lsT, c, newElemTags[2][0], borderElemTags[1]);
int physTag = getBorderTag(lsT, c, newPhysTags[2][0], borderPhysTags[1]);
std::vector<int> phys; phys.push_back(physTag);
@@ -749,10 +864,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
for(int j = 0; j < 3; j++){
int numV = getElementVertexNum(triangles[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j), triangles[i]->z(j), 0);
+ MVertex *newv = new MVertex(triangles[i]->x(j), triangles[i]->y(j), triangles[i]->z(j));
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
- if (!it.second) delete newv;
+ if (!it.second) newv->deleteLast();
}
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
@@ -802,10 +917,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
for(int j = 0; j < 2; j++){
int numV = getElementVertexNum(lines[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j), 0);
+ MVertex *newv = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j));
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
- if (!it.second) delete newv;
+ if (!it.second) newv->deleteLast();
}
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
@@ -822,7 +937,7 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
else lin = new MLine(mv[0], mv[1], ++numEle, ePart);
int lsL = lines[i]->lsTag();
int c = elements[1].count(lsL);
- // suppose that the lines have been cut before the surfaces!
+ // the lines are cut before the surfaces!
int reg = getBorderTag(lsL, c, newElemTags[1][0], borderElemTags[0]);
int physTag = getBorderTag(lsL, c, newPhysTags[1][0], borderPhysTags[0]);
std::vector<int> phys; phys.push_back(physTag);
@@ -844,10 +959,10 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
for(int j = 0; j < 2; j++){
int numV = getElementVertexNum(lines[i]->pt(j), e);
if(numV == -1) {
- MVertex *newv = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j), 0);
+ MVertex *newv = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j));
std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
mv[j] = *(it.first);
- if (!it.second) delete newv;
+ if (!it.second) newv->deleteLast();
}
else {
std::map<int, MVertex*>::iterator it = vertexMap.find(numV);
@@ -903,19 +1018,32 @@ static void elementCutMesh(MElement *e, std::vector<const gLevelset *> &RPN,
GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::map<int, std::vector<MElement*> > elements[10],
std::map<int, MVertex*> &vertexMap,
- std::map<int, std::map<int, std::string> > physicals[4])
+ std::map<int, std::map<int, std::string> > physicals[4],
+ bool cutElem)
{
#if defined(HAVE_DINTEGRATION)
- int numVert = gm->indexMeshVertices(true);
GModel *cutGM = new GModel(gm->getName() + "_cut");
cutGM->setFileName(cutGM->getName());
- newVerticesContainer newVertices;
-
std::vector<GEntity*> gmEntities;
gm->getEntities(gmEntities);
- int numEle = gm->getNumMeshElements();
+
+ std::vector<const gLevelset *> primitives;
+ ls->getPrimitives(primitives);
+ int numVert = gm->indexMeshVertices(true);
+ int nbLs = (cutElem) ? primitives.size() : 1;
+ fullMatrix<double> verticesLs(numVert + 1, nbLs);
+ for(unsigned int i = 0; i < gmEntities.size(); i++) {
+ for(unsigned int j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
+ MVertex *vi = gmEntities[i]->getMeshVertex(j);
+ if(cutElem)
+ for(unsigned int k = 0; k < primitives.size(); k++)
+ verticesLs(vi->getIndex(), k) = (*primitives[k])(vi->x(), vi->y(), vi->z());
+ else
+ verticesLs(vi->getIndex(), 0) = (*ls)(vi->x(), vi->y(), vi->z());
+ }
+ }
std::map<int, int> newElemTags[4];
std::map<int, int> newPhysTags[4];
@@ -923,21 +1051,24 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
newElemTags[d][0] = gm->getMaxElementaryNumber(d);
newPhysTags[d][0] = gm->getMaxPhysicalNumber(d);
}
- std::map<int, int> borderElemTags[2];
- std::map<int, int> borderPhysTags[2];
+ int numEle = gm->getNumMeshElements();
- std::vector<const gLevelset *> primitives;
- ls->getPrimitives(primitives);
- fullMatrix<double> verticesLs(numVert + 1, primitives.size());
- for(unsigned int i = 0; i < gmEntities.size(); i++) {
- for(unsigned int j = 0; j < gmEntities[i]->getNumMeshVertices(); j++) {
- MVertex *vi = gmEntities[i]->getMeshVertex(j);
- for(unsigned int k = 0; k < primitives.size(); k++) {
- verticesLs(vi->getIndex(), k) = (*primitives[k])(vi->x(), vi->y(), vi->z());
+ if(!cutElem) {
+ 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();
+ elementSplitMesh(e, verticesLs, gmEntities[i], gm, numEle, vertexMap,
+ elements, physicals, newElemTags, newPhysTags);
+ cutGM->getMeshPartitions().insert(e->getPartition());
}
}
+ return cutGM;
}
+ newVerticesContainer newVertices;
+ std::map<int, int> borderElemTags[2];
+ std::map<int, int> borderPhysTags[2];
DI_Point::Container cp;
std::vector<DI_Line *> lines;
std::vector<DI_Triangle *> triangles;
diff --git a/Geo/MElementCut.h b/Geo/MElementCut.h
index 7d0d2a82ea..c58dc3202a 100644
--- a/Geo/MElementCut.h
+++ b/Geo/MElementCut.h
@@ -311,6 +311,7 @@ class MLineBorder : public MLine {
GModel *buildCutMesh(GModel *gm, gLevelset *ls,
std::map<int, std::vector<MElement*> > elements[10],
std::map<int, MVertex*> &vertexMap,
- std::map<int, std::map<int, std::string> > physicals[4]);
+ std::map<int, std::map<int, std::string> > physicals[4],
+ bool cutElem);
#endif
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 6ac10df95c..17506c4b54 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -46,6 +46,7 @@ class MVertex{
public:
MVertex(double x, double y, double z, GEntity *ge=0, int num=0);
virtual ~MVertex(){}
+ inline void deleteLast() {if(_num == _globalNum) _globalNum--; delete this;}
// get/reset the global node number
static int getGlobalNumber(){ return _globalNum; }
diff --git a/Numeric/fullMatrix.h b/Numeric/fullMatrix.h
index dcf652c6ca..a159613bff 100644
--- a/Numeric/fullMatrix.h
+++ b/Numeric/fullMatrix.h
@@ -89,10 +89,10 @@ class fullVector
printf("\n");
}
void binarySave (FILE *f) const{
- fwrite (_data,sizeof(scalar),_r,f);
+ fwrite (_data, sizeof(scalar), _r, f);
}
void binaryLoad (FILE *f){
- fread (_data,sizeof(scalar),_r,f);
+ if(fread (_data, sizeof(scalar), _r, f) != _r) return;
}
};
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index a2c0e3a977..28b904abac 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -2009,6 +2009,12 @@ LevelSet :
GM->buildCutGModel(FindLevelSet(t)->ls);
GM->setVisibility(0);
}
+ else if(!strcmp($2, "SplitMesh")){
+ int t = (int)$4;
+ GModel *GM = GModel::current();
+ GM->buildCutGModel(FindLevelSet(t)->ls, false);
+ GM->setVisibility(0);
+ }
else
yymsg(0, "Wrong levelset definition");
Free($2);
diff --git a/contrib/DiscreteIntegration/Integration3D.h b/contrib/DiscreteIntegration/Integration3D.h
index 2665373bb5..f8a8afa14d 100644
--- a/contrib/DiscreteIntegration/Integration3D.h
+++ b/contrib/DiscreteIntegration/Integration3D.h
@@ -43,11 +43,8 @@ class DI_Point
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);}
- // DI_Point(const DI_Point &p) : x_(p.x()), y_(p.y()), z_(p.z()) {Ls.clear(); Ls = p.Ls;}
// assignment
DI_Point & operator=(const DI_Point & rhs);
- // destructor
- // virtual ~DI_Point () {}
// add a levelset value (adjusted to 0 if ls<ZERO_LS_TOL) into the vector Ls
void addLs (const double ls);
// add a levelset value evaluated into the element e
--
GitLab