From f148a0fe10985d9f0b6e3803d9b2b8f1c0e250ea Mon Sep 17 00:00:00 2001
From: Gaetan Bricteux <gaetan.bricteux@uclouvain.be>
Date: Tue, 18 Jun 2013 08:48:38 +0000
Subject: [PATCH] makeSimplyConnected
---
Geo/GModel.cpp | 169 ++++++++++++++++++++++++--------------------
Geo/MElementCut.cpp | 72 ++++++++++++-------
2 files changed, 140 insertions(+), 101 deletions(-)
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 63f4924cd9..b75edfc705 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2304,42 +2304,52 @@ void makeSimplyConnected(std::map<int, std::vector<MElement*> > elements[10])
for(unsigned int j = 0; j < elements[4][ri].size(); j++)
allElements.push_back(elements[4][ri][j]);
std::vector<std::vector<MElement*> > conRegions;
- int nbRegions = connectedVolumes(allElements, conRegions);
- printf("%d regions (%d)\n",nbRegions,ri);
- for(int j = 0; j < nbRegions; j++){
+ int nbConRegions = connectedVolumes(allElements, conRegions);
+ for(int j = 0; j < nbConRegions; j++){
if(conRegions[j].size() < 3){ //remove conRegions containing 1 or 2 elements
+ //find adjacent region
+ int r2 = ri;
+ if(regs.size() == 2)
+ r2 = (ri + 1) % 2;
+ else{
+ for(unsigned int k = 0; k < conRegions[j].size(); k++){
+ MElement *el = conRegions[j][k];
+ for(int l = 0; l < el->getNumFaces(); l++){
+ MFace mf = el->getFace(l);
+ std::multimap<MFace, MElement*, Less_Face>::iterator itl = f2e.lower_bound(mf);
+ for(; itl != f2e.upper_bound(mf); itl++){
+ if(itl->second != el) break;
+ }
+ MElement *el2 = itl->second;
+ bool sameRegion = false;
+ for(unsigned int m = 0; m < conRegions[j].size(); m++)
+ if(conRegions[j][m] == el2) {
+ sameRegion = true; break;
+ }
+ if(sameRegion) continue;
+ for(unsigned int m = 0; m < regs.size(); m++){
+ int rm = regs[m];
+ if(rm == ri) continue;
+ for(unsigned int n = 0; n < elements[4][rm].size(); n++)
+ if(elements[4][rm][n] == el2){
+ r2 = rm;
+ break;
+ }
+ if(r2 != ri) break;
+ }
+ if(r2 != ri) break;
+ }
+ if(r2 != ri) break;
+ }
+ if(r2 == ri) Msg::Warning("Element not found for simply connected regions");
+ }
+
for(unsigned int k = 0; k < conRegions[j].size(); k++){
MElement *el = conRegions[j][k];
unsigned int l = 0;
for(; l < elements[4][ri].size(); l++)
if(elements[4][ri][l] == el) break;
elements[4][ri].erase(elements[4][ri].begin() + l);
- //find adjacent region
- int r2 = ri;
- if(regs.size() == 2)
- r2 = (ri + 1) % 2;
- else{
- MFace mf = el->getFace(0);
- std::multimap<MFace, MElement*, Less_Face>::iterator itl = f2e.lower_bound(mf);
- for(; itl != f2e.upper_bound(mf); itl++){
- if(itl->second != el) break;
- }
- if(itl == f2e.upper_bound(mf)) printf("adjacent 3d element not found\n");
- MElement *el2 = itl->second;
- bool found = false;
- for(unsigned int m = 0; m < regs.size(); m++){
- int rm = regs[m];
- if(rm == ri) continue;
- for(unsigned int n = 0; n < elements[4][rm].size(); n++)
- if(elements[4][rm][n] == el2){
- r2 = rm;
- found = true;
- break;
- }
- if(found) break;
- }
- if(r2 == ri) Msg::Warning("Element not found for simply connected regions");
- }
elements[4][r2].push_back(el);
}
}
@@ -2367,41 +2377,50 @@ void makeSimplyConnected(std::map<int, std::vector<MElement*> > elements[10])
allElements.push_back(elements[2][fi][j]);
std::vector<std::vector<MElement*> > conSurfaces;
int nbSurfaces = connectedSurfaces(allElements, conSurfaces);
- printf("%d consurfaces (reg=%d)\n",nbSurfaces,fi);
for(int j = 0; j < nbSurfaces; j++){
if(conSurfaces[j].size() < 3){ //remove conSurfaces containing 1 or 2 elements
+ //find adjacent surface
+ int f2 = fi;
+ if(faces.size() == 2)
+ f2 = (fi + 1) % 2;
+ else{
+ for(unsigned int k = 0; k < conSurfaces[j].size(); k++){
+ MElement *el = conSurfaces[j][k];
+ for(int l = 0; l < el->getNumEdges(); l++){
+ MEdge me = el->getEdge(l);
+ std::multimap<MEdge, MElement*, Less_Edge>::iterator itl = e2e.lower_bound(me);
+ for(; itl != e2e.upper_bound(me); itl++){
+ if(itl->second != el) break;
+ }
+ MElement *el2 = itl->second;
+ bool sameSurface = false;
+ for(unsigned int m = 0; m < conSurfaces[j].size(); m++)
+ if(conSurfaces[j][m] == el2) {
+ sameSurface = true; break;
+ }
+ if(sameSurface) continue;
+ for(unsigned int m = 0; m < faces.size(); m++){
+ int fm = faces[m];
+ if(fm == fi) continue;
+ for(unsigned int n = 0; n < elements[2][fm].size(); n++)
+ if(elements[2][fm][n] == el2){
+ f2 = fm;
+ break;
+ }
+ if(f2 != fi) break;
+ }
+ if(f2 != fi) break;
+ }
+ if(f2 != fi) break;
+ }
+ if(f2 == fi) Msg::Warning("Element not found for simply connected surfaces");
+ }
for(unsigned int k = 0; k < conSurfaces[j].size(); k++){
MElement *el = conSurfaces[j][k];
unsigned int l = 0;
for(; l < elements[2][fi].size(); l++)
if(elements[2][fi][l] == el) break;
elements[2][fi].erase(elements[2][fi].begin() + l);
- //find adjacent surface
- int f2 = fi;
- if(faces.size() == 2)
- f2 = (fi + 1) % 2;
- else{
- MEdge me = el->getEdge(0);
- std::multimap<MEdge, MElement*, Less_Edge>::iterator itl = e2e.lower_bound(me);
- for(; itl != e2e.upper_bound(me); itl++){
- if(itl->second != el) break;
- }
- if(itl == e2e.upper_bound(me)) printf("adjacent 2d element not found\n");
- MElement *el2 = itl->second;
- bool found = false;
- for(unsigned int m = 0; m < faces.size(); m++){
- int fm = faces[m];
- if(fm == fi) continue;
- for(unsigned int n = 0; n < elements[2][fm].size(); n++)
- if(elements[2][fm][n] == el2){
- f2 = fm;
- found = true;
- break;
- }
- if(found) break;
- }
- if(f2 == fi) Msg::Warning("Element not found for simply connected surfaces");
- }
elements[2][f2].push_back(el);
}
}
@@ -2427,7 +2446,8 @@ GModel *GModel::buildCutGModel(gLevelset *ls, bool cutElem, bool saveTri)
GModel *cutGM = buildCutMesh(this, ls, elements, vertexMap, physicals, cutElem);
- makeSimplyConnected(elements);
+ if(!cutElem)
+ makeSimplyConnected(elements);
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
cutGM->_storeElementsInEntities(elements[i]);
@@ -2670,15 +2690,15 @@ GRegion* GModel::addVolume (std::vector<std::vector<GFace *> > faces){
GFace *GModel::add2Drect(double x0, double y0, double dx, double dy)
{
- if(_factory)
- return _factory->add2Drect(this, x0, y0, dx, dy);
+ if(_factory)
+ return _factory->add2Drect(this, x0, y0, dx, dy);
return 0;
}
GFace *GModel::add2Dellips(double xc, double yc, double rx, double ry)
{
- if(_factory)
- return _factory->add2Dellips(this, xc, yc, rx, ry);
+ if(_factory)
+ return _factory->add2Dellips(this, xc, yc, rx, ry);
return 0;
}
@@ -2740,8 +2760,7 @@ GEntity *GModel::addBlock(std::vector<double> p1, std::vector<double> p2)
GEntity *GModel::add3DBlock(std::vector<double> p1, double dx, double dy, double dz )
{
- if(_factory)
- return _factory->add3DBlock(this, p1, dx, dy, dz);
+ if(_factory) return _factory->add3DBlock(this, p1, dx, dy, dz);
return 0;
}
@@ -2775,32 +2794,32 @@ GModel *GModel::computeBooleanDifference(GModel *tool, int createNewModel)
void GModel::salomeconnect()
{
- if(_factory)
- _factory->salomeconnect(this);
+ if(_factory) _factory->salomeconnect(this);
}
void GModel::occconnect()
{
- if(_factory)
- _factory->occconnect(this);
+ if(_factory) _factory->occconnect(this);
}
void GModel::setPeriodicAllFaces(std::vector<double> FaceTranslationVector)
{
- if(_factory)
- _factory->setPeriodicAllFaces(this, FaceTranslationVector);
+ if(_factory) _factory->setPeriodicAllFaces(this, FaceTranslationVector);
}
-void GModel::setPeriodicPairOfFaces(int numFaceMaster, std::vector<int> EdgeListMaster, int numFaceSlave, std::vector<int> EdgeListSlave)
+void GModel::setPeriodicPairOfFaces(int numFaceMaster, std::vector<int> EdgeListMaster,
+ int numFaceSlave, std::vector<int> EdgeListSlave)
{
- if(_factory)
- _factory->setPeriodicPairOfFaces(this, numFaceMaster, EdgeListMaster, numFaceSlave, EdgeListSlave);
+ if(_factory)
+ _factory->setPeriodicPairOfFaces(this, numFaceMaster, EdgeListMaster,
+ numFaceSlave, EdgeListSlave);
}
-void GModel::setPhysicalNumToEntitiesInBox(int EntityType, int PhysicalGroupNumber, std::vector<double> p1, std::vector<double> p2)
+void GModel::setPhysicalNumToEntitiesInBox(int EntityType, int PhysicalGroupNumber,
+ std::vector<double> p1, std::vector<double> p2)
{
- if(_factory)
- _factory->setPhysicalNumToEntitiesInBox(this,EntityType,PhysicalGroupNumber,p1,p2);
+ if(_factory)
+ _factory->setPhysicalNumToEntitiesInBox(this, EntityType, PhysicalGroupNumber, p1, p2);
}
static void computeDuplicates(GModel *model,
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index 44a5808738..744646ee64 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -557,7 +557,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::map<int, int> &newPhysTags)
+static void getPhysicalTag(int lsTag, const std::vector<int> &physicals,
+ std::map<int, int> &newPhysTags)
{
for(unsigned int i = 0; i < physicals.size(); i++){
int phys = physicals[i];
@@ -605,7 +606,7 @@ static void elementSplitMesh(MElement *e, std::vector<gLevelset *> &RPN,
int gLsTag = RPN[RPN.size() - 1]->getTag();
MElement *copy = e->copy(vertexMap, newParents, newDomains);
- double eps = 0.0; //-1.e-6;
+ double eps = 1.e-6;
bool splitElem = false;
//split acording to center of gravity
@@ -766,7 +767,8 @@ static void elementSplitMesh(MElement *e, std::vector<gLevelset *> &RPN,
if(f1 == f2) break;
}
if(i < 0){
- MTriangle *tri = new MTriangle(f1.getVertex(0), f1.getVertex(1), f1.getVertex(2));
+ MTriangle *tri = new MTriangle(f1.getVertex(0), f1.getVertex(1),
+ f1.getVertex(2));
elements[2][reg].push_back(tri);
}
else{
@@ -923,7 +925,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
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());
+ 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, recur, nodeLs);
isCut = iC1 || iC2 || iC3;
@@ -944,7 +947,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
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());
+ 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, recur, nodeLs);
isCut = iC1 || iC2;
@@ -973,8 +977,10 @@ static void elementCutMesh(MElement *e, std::vector<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));
- std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
+ 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) newv->deleteLast();
}
@@ -1013,7 +1019,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
}
// check for border surfaces cut earlier along the polyhedra
std::pair<std::multimap<MElement*, MElement*>::iterator,
- std::multimap<MElement*, MElement*>::iterator> itr = borders[1].equal_range(copy);
+ std::multimap<MElement*, MElement*>::iterator> itr =
+ borders[1].equal_range(copy);
std::vector<std::pair<MElement*, MElement*> > bords;
for(std::multimap<MElement*, MElement*>::iterator it = itr.first;
it != itr.second; it++) {
@@ -1060,7 +1067,8 @@ static void elementCutMesh(MElement *e, std::vector<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));
+ 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) newv->deleteLast();
@@ -1082,7 +1090,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
}
else tri = new MTriangle(mv[0], mv[1], mv[2], ++numEle, ePart);
int lsTag = triangles[i]->lsTag();
- int cR = elements[2].count(lsTag) + elements[3].count(lsTag) + elements[8].count(lsTag);
+ int cR = elements[2].count(lsTag) + elements[3].count(lsTag) +
+ elements[8].count(lsTag);
int cP = 0;
for(std::map<int, std::map<int, std::string> >::iterator it = physicals[2].begin();
it != physicals[2].end(); it++)
@@ -1152,8 +1161,10 @@ static void elementCutMesh(MElement *e, std::vector<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));
- std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
+ 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) newv->deleteLast();
}
@@ -1180,7 +1191,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
int numV = getElementVertexNum(quads[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);
+ std::pair<newVerticesContainer::iterator, bool> it =
+ newVertices.insert(newv);
mv[j] = *(it.first);
if (!it.second) newv->deleteLast();
}
@@ -1237,7 +1249,8 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
}
// check for border lines cut earlier along the polygons
std::pair<std::multimap<MElement*, MElement*>::iterator,
- std::multimap<MElement*, MElement*>::iterator> itr = borders[0].equal_range(copy);
+ std::multimap<MElement*, MElement*>::iterator> itr =
+ borders[0].equal_range(copy);
std::vector<std::pair<MElement*, MElement*> > bords;
for(std::multimap<MElement*, MElement*>::iterator it = itr.first;
it != itr.second; ++it) {
@@ -1342,14 +1355,16 @@ static void elementCutMesh(MElement *e, std::vector<gLevelset *> &RPN,
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));
- std::pair<newVerticesContainer::iterator, bool> it = newVertices.insert(newv);
+ 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(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j), 0, numV);
+ mv[j] = new MVertex(lines[i]->x(j), lines[i]->y(j), lines[i]->z(j),
+ 0, numV);
vertexMap[numV] = mv[j];
}
else mv[j] = it->second;
@@ -1458,7 +1473,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
}
}
- int numEle = gm->getNumMeshElements() + gm->getNumMeshParentElements(); //element number increment
+ //element number increment
+ int numEle = gm->getNumMeshElements() + gm->getNumMeshParentElements();
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);
@@ -1488,9 +1504,9 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
for(unsigned int j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
MElement *e = gmEntities[i]->getMeshElement(j);
e->setVolumePositive();
- elementSplitMesh(e, RPN, verticesLs, gmEntities[i], gm, numEle, vertexMap, newParents,
- newDomains, elements, physicals, newElemTags, newPhysTags,
- borderElemTags, borderPhysTags);
+ elementSplitMesh(e, RPN, verticesLs, gmEntities[i], gm, numEle, vertexMap,
+ newParents, newDomains, elements, physicals, newElemTags,
+ newPhysTags, borderElemTags, borderPhysTags);
cutGM->getMeshPartitions().insert(e->getPartition());
}
}
@@ -1499,7 +1515,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
//CutMesh
newVerticesContainer newVertices;
- std::multimap<MElement*, MElement*> borders[2]; //multimap<domain,border>[polyg=0,polyh=1]
+ //multimap<domain,border>[polyg=0,polyh=1]
+ std::multimap<MElement*, MElement*> borders[2];
DI_Point::Container cp;
std::vector<DI_Line *> lines;
std::vector<DI_Triangle *> triangles;
@@ -1516,9 +1533,10 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
for(unsigned int j = 0; j < gmEntities[i]->getNumMeshElements(); j++) {
MElement *e = gmEntities[i]->getMeshElement(j);
e->setVolumePositive();
- elementCutMesh(e, RPN, verticesLs, gmEntities[i], gm, numEle, vertexMap, newVertices, newParents,
- newDomains, borders, elements, physicals, newElemTags, newPhysTags,
- borderElemTags, borderPhysTags, cp, lines, triangles, quads, tetras, hexas);
+ elementCutMesh(e, RPN, verticesLs, gmEntities[i], gm, numEle, vertexMap,
+ newVertices, newParents, newDomains, borders, elements, physicals,
+ newElemTags, newPhysTags, borderElemTags, borderPhysTags, cp,
+ lines, triangles, quads, tetras, hexas);
cutGM->getMeshPartitions().insert(e->getPartition());
}
@@ -1581,7 +1599,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
for(unsigned int k = 0; k < quads.size(); k++) delete quads[k];
for(unsigned int k = 0; k < tetras.size(); k++) delete tetras[k];
for(unsigned int k = 0; k < hexas.size(); k++) delete hexas[k];
- cp.clear(); lines.clear(); triangles.clear(); quads.clear(); tetras.clear(); hexas.clear();
+ cp.clear(); lines.clear(); triangles.clear(); quads.clear();
+ tetras.clear(); hexas.clear();
}
#if 0
@@ -1608,7 +1627,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
printf("\n");
#endif
- for(newVerticesContainer::iterator it = newVertices.begin() ; it != newVertices.end(); ++it) {
+ for(newVerticesContainer::iterator it = newVertices.begin();
+ it != newVertices.end(); ++it){
vertexMap[(*it)->getNum()] = *it;
}
--
GitLab