From 0b3bfed5472c09536662827b8d33a2eb5f273953 Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Wed, 19 May 2010 15:31:35 +0000
Subject: [PATCH] Spectral conformal faster + createTopo for connected faces
done !
---
Geo/GFaceCompound.cpp | 45 +++------
Geo/GModel.cpp | 185 ++++++++++++++++++++++++-----------
Geo/GModelIO_Mesh.cpp | 4 +-
Mesh/multiscalePartition.cpp | 2 +-
Solver/multiscaleLaplace.cpp | 1 -
benchmarks/boolean/nurbs.lua | 2 +-
6 files changed, 147 insertions(+), 92 deletions(-)
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 460d0099b0..1c4ae412ba 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -363,7 +363,7 @@ bool GFaceCompound::trivial() const
//of triangles
bool GFaceCompound::checkFolding(std::vector<MVertex*> &ordered) const
{
-
+
bool has_no_folding = true;
for(unsigned int i = 0; i < ordered.size()-1; ++i){
@@ -519,10 +519,10 @@ bool GFaceCompound::parametrize() const
else if (_mapping == CONFORMAL){
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- //bool withoutFolding = parametrize_conformal_spectral() ;
- bool withoutFolding = parametrize_conformal();
+ bool withoutFolding = parametrize_conformal_spectral() ;
+ //bool withoutFolding = parametrize_conformal();
if ( withoutFolding == false ){
- //printStuff(); exit(1);
+ printStuff(); //exit(1);
Msg::Warning("$$$ Parametrization switched to harmonic map");
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
@@ -1129,6 +1129,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
linearSystem <double> *lsysB = new linearSystemPETSc<double>;
dofManager<double> myAssembler(lsysA, lsysB);
+ //printf("nbNodes = %d in spectral param \n", allNodes.size());
//-------------------------------
myAssembler.setCurrentMatrix("A");
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -1164,26 +1165,6 @@ bool GFaceCompound::parametrize_conformal_spectral() const
//-------------------------------
myAssembler.setCurrentMatrix("B");
-
-// massTerm mass1(model(), 1, &MONE);
-// massTerm mass2(model(), 2, &MONE);
-// it = _compound.begin();
-// for( ; it != _compound.end() ; ++it){
-// for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
-// SElement se((*it)->triangles[i]);
-// mass1.addToMatrix(myAssembler, &se);
-// mass2.addToMatrix(myAssembler, &se);
-// }
-// }
-
-// std::list<GEdge*>::const_iterator it0 = _U0.begin();
-// for( ; it0 != _U0.end(); ++it0 ){
-// for(unsigned int i = 0; i < (*it0)->lines.size(); i++ ){
-// SElement se((*it0)->lines[i]);
-// mass1.addToMatrix(myAssembler, &se);
-// mass2.addToMatrix(myAssembler, &se);
-// }
-// }
double small = 0.0;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
@@ -1198,19 +1179,25 @@ bool GFaceCompound::parametrize_conformal_spectral() const
}
}
+ //mettre max NC contraintes par bord
int NB = ordered.size();
- for (int i = 0; i< NB; i++){
- MVertex *v1 = ordered[i];
- for (int j = 0; j< NB; j++){
- MVertex *v2 = ordered[j];
+ int NC = std::min(100,NB);
+ int jump = (int) NB/NC;
+ int nbLoop = (int) NB/jump ;
+ //printf("nb bound nodes=%d jump =%d \n", NB, jump);
+ for (int i = 0; i< nbLoop; i++){
+ MVertex *v1 = ordered[i*jump];
+ for (int j = 0; j< nbLoop; j++){
+ MVertex *v2 = ordered[j*jump];
myAssembler.assemble(v1, 0, 1, v2, 0, 1, -1./NB);
myAssembler.assemble(v1, 0, 2, v2, 0, 2, -1./NB);
}
}
//-------------------------------
+ //printf("Solve eigensystem \n");
eigenSolver eig(&myAssembler, "B" , "A", true);
- int nb = 1;
+ int nb = 2;
bool converged = eig.solve(nb, "largest");
if(converged) {
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 01f8e88eb8..d2997b089a 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -40,6 +40,90 @@
std::vector<GModel*> GModel::list;
int GModel::_current = -1;
+static void recur_connect(MVertex *v,
+ std::multimap<MVertex*,MEdge> &v2e,
+ std::set<MEdge,Less_Edge> &group,
+ std::set<MVertex*> &touched)
+{
+ if (touched.find(v) != touched.end()) return;
+
+ touched.insert(v);
+ for (std::multimap <MVertex*,MEdge>::iterator it = v2e.lower_bound(v);
+ it != v2e.upper_bound(v) ; ++it){
+ group.insert(it->second);
+ for (int i = 0; i < it->second.getNumVertices(); ++i){
+ recur_connect (it->second.getVertex(i), v2e, group, touched);
+ }
+ }
+}
+
+// starting form a list of elements, returns lists of lists that are
+// all simply connected
+static void recur_connect_e(const MEdge &e,
+ std::multimap<MEdge,MElement*, Less_Edge> &e2e,
+ std::set<MElement*> &group,
+ std::set<MEdge, Less_Edge> &touched)
+{
+ if (touched.find(e) != touched.end())return;
+ touched.insert(e);
+ for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.lower_bound(e);
+ it != e2e.upper_bound(e) ; ++it){
+ group.insert(it->second);
+ for (int i = 0; i < it->second->getNumEdges(); ++i){
+ recur_connect_e(it->second->getEdge(i), e2e, group, touched);
+ }
+ }
+}
+
+static int connected_bounds(std::set<MEdge, Less_Edge> &edges,
+ std::vector<std::vector<MEdge> > &boundaries)
+{
+ std::multimap<MVertex*,MEdge> v2e;
+ for(std::set<MEdge, Less_Edge>::iterator it = edges.begin(); it != edges.end(); it++){
+ for (int j=0;j<it->getNumVertices();j++){
+ v2e.insert(std::make_pair(it->getVertex(j),*it));
+ }
+ }
+
+ while (!v2e.empty()){
+ std::set<MEdge, Less_Edge> group;
+ std::set<MVertex*> touched;
+ recur_connect(v2e.begin()->first, v2e, group, touched);
+ std::vector<MEdge> temp;
+ temp.insert(temp.begin(), group.begin(), group.end());
+ boundaries.push_back(temp);
+ for (std::set<MVertex*>::iterator it = touched.begin() ; it != touched.end();++it)
+ v2e.erase(*it);
+ }
+
+ return boundaries.size();
+}
+
+//--------------------------------------------------------------
+static int connectedRegions (std::vector<MElement*> &elements,
+ std::vector<std::vector<MElement*> > ®ions)
+{
+ std::multimap<MEdge,MElement*,Less_Edge> e2e;
+ for (unsigned int i = 0; i < elements.size(); ++i){
+ for (int j = 0; j < elements[i]->getNumEdges(); j++){
+ e2e.insert(std::make_pair(elements[i]->getEdge(j),elements[i]));
+ }
+ }
+ while (!e2e.empty()){
+ std::set<MElement*> group;
+ std::set<MEdge,Less_Edge> touched;
+ recur_connect_e (e2e.begin()->first,e2e,group,touched);
+ std::vector<MElement*> temp;
+ temp.insert(temp.begin(), group.begin(), group.end());
+ regions.push_back(temp);
+ for ( std::set<MEdge,Less_Edge>::iterator it = touched.begin() ; it != touched.end();++it)
+ e2e.erase(*it);
+ }
+
+ return regions.size();
+}
+//-----------------------------------------------------------------------------
+
GModel::GModel(std::string name)
: _name(name), _visible(1), _octree(0),
_geo_internals(0), _occ_internals(0), _acis_internals(0), _fm_internals(0),_factory(0), _fields(0), _currentMeshEntity(0), normals(0)
@@ -1079,7 +1163,7 @@ void GModel::createTopologyFromMesh()
discRegions.push_back((discreteRegion*) *it);
//EMI-FIX in case of createTopology for Volumes
- //all faces are set to aeach volume
+ //all faces are set to each volume
for (std::vector<discreteRegion*>::iterator it = discRegions.begin();
it != discRegions.end(); it++)
(*it)->setBoundFaces();
@@ -1090,8 +1174,48 @@ void GModel::createTopologyFromMesh()
if((*it)->geomType() == GEntity::DiscreteSurface)
discFaces.push_back((discreteFace*) *it);
+ //--------
//EMI TODO
//check for closed faces
+
+ printf("discr faces size=%d \n", discFaces.size());
+ std::vector<discreteFace*> newDiscFaces;
+ for (std::vector<discreteFace*>::iterator itF = discFaces.begin();
+ itF != discFaces.end(); itF++){
+
+ std::vector<MElement*> allElements;
+ for (unsigned int i = 0; i < (*itF)->getNumMeshElements(); i++)
+ allElements.push_back((*itF)->getMeshElement(i));
+
+ std::vector<std::vector<MElement*> > conRegions;
+ int nbFaces = connectedRegions (allElements, conRegions);
+
+ remove(*itF);
+ for (int ifa = 0; ifa < nbFaces; ifa++){
+ std::vector<MElement*> myElements = conRegions[ifa];
+
+ int numF = maxFaceNum()+1;
+ discreteFace *f = new discreteFace(this, numF);
+ //printf("*** Created discreteFace %d \n", numF);
+ add(f);
+ newDiscFaces.push_back(f);
+
+ //fill new face with mesh vertices
+ std::set<MVertex*> allV;
+ for(unsigned int i = 0; i < myElements.size(); i++) {
+ MVertex *v0 = myElements[i]->getVertex(0);
+ MVertex *v1 = myElements[i]->getVertex(1);
+ MVertex *v2 = myElements[i]->getVertex(2);
+ f->triangles.push_back(new MTriangle( v0, v1, v2));
+ allV.insert(v0); allV.insert(v1); allV.insert(v2);
+ v0->setEntity(f); v1->setEntity(f); v2->setEntity(f);
+ }
+ f->mesh_vertices.insert(f->mesh_vertices.begin(), allV.begin(),allV.end());
+ }
+
+ }
+ discFaces = newDiscFaces;
+ //------
createTopologyFromFaces(discFaces);
@@ -1102,64 +1226,7 @@ void GModel::createTopologyFromMesh()
exportDiscreteGEOInternals();
}
-static void recur_connect(MVertex *v,
- std::multimap<MVertex*,MEdge> &v2e,
- std::set<MEdge,Less_Edge> &group,
- std::set<MVertex*> &touched)
-{
- if (touched.find(v) != touched.end()) return;
-
- touched.insert(v);
- for (std::multimap <MVertex*,MEdge>::iterator it = v2e.lower_bound(v);
- it != v2e.upper_bound(v) ; ++it){
- group.insert(it->second);
- for (int i = 0; i < it->second.getNumVertices(); ++i){
- recur_connect (it->second.getVertex(i), v2e, group, touched);
- }
- }
-}
-
-// starting form a list of elements, returns lists of lists that are
-// all simply connected
-static void recur_connect_e(const MEdge &e,
- std::multimap<MEdge,MElement*, Less_Edge> &e2e,
- std::set<MElement*> &group,
- std::set<MEdge, Less_Edge> &touched)
-{
- if (touched.find(e) != touched.end())return;
- touched.insert(e);
- for (std::multimap <MEdge,MElement*,Less_Edge>::iterator it = e2e.lower_bound(e);
- it != e2e.upper_bound(e) ; ++it){
- group.insert(it->second);
- for (int i = 0; i < it->second->getNumEdges(); ++i){
- recur_connect_e(it->second->getEdge(i), e2e, group, touched);
- }
- }
-}
-
-static int connected_bounds(std::set<MEdge, Less_Edge> &edges,
- std::vector<std::vector<MEdge> > &boundaries)
-{
- std::multimap<MVertex*,MEdge> v2e;
- for(std::set<MEdge, Less_Edge>::iterator it = edges.begin(); it != edges.end(); it++){
- for (int j=0;j<it->getNumVertices();j++){
- v2e.insert(std::make_pair(it->getVertex(j),*it));
- }
- }
-
- while (!v2e.empty()){
- std::set<MEdge, Less_Edge> group;
- std::set<MVertex*> touched;
- recur_connect(v2e.begin()->first, v2e, group, touched);
- std::vector<MEdge> temp;
- temp.insert(temp.begin(), group.begin(), group.end());
- boundaries.push_back(temp);
- for (std::set<MVertex*>::iterator it = touched.begin() ; it != touched.end();++it)
- v2e.erase(*it);
- }
- return boundaries.size();
-}
void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
{
@@ -1243,7 +1310,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
}
std::vector<std::vector<MEdge> > boundaries;
- double nbBounds = connected_bounds(myEdges, boundaries);
+ int nbBounds = connected_bounds(myEdges, boundaries);
for (int ib = 0; ib < nbBounds; ib++){
std::vector<MEdge> myLines = boundaries[ib];
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index f6771098d1..f91c9873e2 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -535,7 +535,9 @@ static void writeElementMSH(FILE *fp, GModel *model, T *ele, bool saveAll,
for(unsigned int j = 0; j < physicals.size(); j++)
ele->writeMSH(fp, version, binary, ++num, elementary, physicals[j],
parentNum, &ghosts);
- (*newElemNumbers)(ele->getNum()) = num;
+
+ //(*newElemNumbers)(ele->getNum()) = num;
+
}
template<class T>
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 2da3da6241..25003d7098 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -348,7 +348,7 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
int genus, AR, NB;
getGenusAndRatio(regions[i], genus, AR, NB);
- printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
+ //printLevel (nextLevel->elements, nextLevel->recur,nextLevel->region);
if (genus < 0) {
Msg::Error("Genus partition is negative G=%d!", genus);
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index a0df0d798d..d23372a6b7 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -1006,7 +1006,6 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
else tooSmall.insert(tooSmall.begin(), regGoodSize[i].begin(), regGoodSize[i].end());
}
}
-
//Add the not too small regions to the level.elements
std::vector<std::vector<MElement*> > regions_, regions ;
diff --git a/benchmarks/boolean/nurbs.lua b/benchmarks/boolean/nurbs.lua
index ef04763fd8..97c3efa30d 100644
--- a/benchmarks/boolean/nurbs.lua
+++ b/benchmarks/boolean/nurbs.lua
@@ -21,7 +21,7 @@ g:addLine(v7,v2);
g:addBezier (v1,v2, {{v3:x(),v3:y(),0},{v4:x(),v4:y(),0},{v5:x(),v5:y(),0},{v6:x(),v6:y(),0},{v7:x(),v7:y(),0}});
g:addNURBS (v1,v2, {{v3:x(),v3:y(),0},{v4:x(),v4:y(),0},{v5:x(),v5:y(),0},{v6:x(),v6:y(),0},{v7:x(),v7:y(),0}},
- {0,0.25,0.5,0.75,1},{1,1,1,1,1,1,1},{4,1,1,1,4});
+ {0,0.25,0.5,0.75,1}}},{1,1,1,1,1,1,1},{4,1,1,1,4});
g:addNURBS (v1,v2, {{v3:x(),v3:y(),0},{v4:x(),v4:y(),0},{v5:x(),v5:y(),0},{v6:x(),v6:y(),0},{v7:x(),v7:y(),0}},
{0,0.5,1},{1,1,1,1,1,1,1},{4,3,4});
--
GitLab