diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index 80360803c32dd200264662ec46b08c853318a10e..e20d694399bfe032302fdb6eb643627df74636b0 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -520,7 +520,7 @@ int MergeFile(const std::string &fileName, bool warnIfMissing, bool setBoundingB
status = GModel::readGEO(fileName);
}
}
-
+
ComputeMaxEntityNum();
if(setBoundingBox) SetBoundingBox();
CTX::instance()->geom.draw = 1;
diff --git a/Geo/CMakeLists.txt b/Geo/CMakeLists.txt
index ecfd656e588526d43ad5c4427716e551268d700c..1c2dc82e20a76cb7c647264138c52cc536bcc8ea 100644
--- a/Geo/CMakeLists.txt
+++ b/Geo/CMakeLists.txt
@@ -21,6 +21,7 @@ set(SRC
ACISEdge.cpp
ACISFace.cpp
GModel.cpp
+ GModelCreateTopologyFromMesh.cpp
GModelVertexArrays.cpp
GModelFactory.cpp
GModelIO_GEO.cpp GModelIO_ACIS.cpp GModelIO_OCC.cpp GModelIO_TOCHNOG.cpp
diff --git a/Geo/GEdge.h b/Geo/GEdge.h
index e0db7f1cb004995a780aef3536a5b7fbe34b598c..d88bbe461a9ca235b9a1c3e723fe985cdb460de3 100644
--- a/Geo/GEdge.h
+++ b/Geo/GEdge.h
@@ -48,6 +48,8 @@ class GEdge : public GEntity{
virtual void deleteMesh();
// get the start/end vertices of the edge
+ void setBeginVertex(GVertex *gv) { v0=gv; }
+ void setEndVertex(GVertex *gv) { v1=gv; }
virtual GVertex *getBeginVertex() const { return v0; }
virtual GVertex *getEndVertex() const { return v1; }
diff --git a/Geo/GFace.h b/Geo/GFace.h
index f394f44fbe347a9b7fdfaef9d89616484a29ac74..b0a89bd11f6bd90666338396b99ba70faf4c321e 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -112,6 +112,7 @@ class GFace : public GEntity{
// edges that bound the face
virtual std::list<GEdge*> edges() const { return l_edges; }
+ inline void set(const std::list<GEdge*> f) { l_edges= f; }
virtual std::list<int> edgeOrientations() const { return l_dirs; }
inline bool containsEdge (int iEdge) const
{
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 4d61c62ccc84f3e14224c1207f0109d7cd841543..1aac4eb9a25f1f9d50325fc04e7e83dd2c7249b7 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -42,6 +42,7 @@
#include "meshGEdge.h"
#include "meshGFace.h"
#include "meshGRegion.h"
+#include "GModelCreateTopologyFromMesh.h"
#if defined(HAVE_MESH)
#include "Field.h"
@@ -2251,15 +2252,26 @@ void GModel::makeDiscreteFacesSimplyConnected()
Msg::Debug("Done making discrete faces simply connected");
}
+
void GModel::createTopologyFromMesh(int ignoreHoles)
{
Msg::StatusBar(true, "Creating topology from mesh...");
double t1 = Cpu();
-
removeDuplicateMeshVertices(CTX::instance()->geom.tolerance);
makeDiscreteRegionsSimplyConnected();
makeDiscreteFacesSimplyConnected();
+ // // TEST !!!!!!!!
+ // if (1)
+ // {
+ // createTopologyFromMeshNew (this);
+ // double t2 = Cpu();
+ // Msg::StatusBar(true, "Done creating topology from mesh (%g s)", t2 - t1);
+ // return;
+ // }
+
+
+
// create topology for all discrete regions
std::vector<discreteRegion*> discRegions;
for(riter it = firstRegion(); it != lastRegion(); it++)
diff --git a/Geo/GModelCreateTopologyFromMesh.cpp b/Geo/GModelCreateTopologyFromMesh.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..24d2cfbe9d5c8e7686c33a0f2680d9f99c591528
--- /dev/null
+++ b/Geo/GModelCreateTopologyFromMesh.cpp
@@ -0,0 +1,328 @@
+#include <map>
+#include "GModelCreateTopologyFromMesh.h"
+#include "GModel.h"
+#include "Geo.h"
+#include "discreteFace.h"
+#include "discreteEdge.h"
+#include "MPoint.h"
+#include "MVertex.h"
+#include "MLine.h"
+#include "MEdge.h"
+#include "MFace.h"
+#include "MTriangle.h"
+#include "MQuadrangle.h"
+
+/*
+ Assumptions :
+ --> The input mesh is potentially (partially) coloured
+ -->
+*/
+
+template <class T>
+class myBundle {
+public :
+ std::set<T> stuff;
+ void insert (T t) {stuff.insert(t);}
+ void print() const {
+ printf("bundle of %ld entities\n",stuff.size());
+ for (typename std::set<T>::iterator it = stuff.begin(); it != stuff.end() ; ++it){
+ printf ("%d ",(*it)->tag());
+ }
+ printf("\n");
+ }
+ bool operator < (const myBundle<T> & other) const{
+ if (other.stuff.size() < stuff.size()) return true;
+ if (other.stuff.size() > stuff.size()) return false;
+ typename std::set<T>::const_iterator it = stuff.begin();
+ typename std::set<T>::const_iterator it2 = other.stuff.begin();
+ for ( ; it != stuff.end() ; ++it, ++it2){
+ if (*it < *it2) return true;
+ if (*it > *it2) return false;
+ }
+ return false;
+ }
+};
+
+
+void createTopologyFromMesh1D ( GModel *gm ) {
+ std::map<MVertex*, myBundle<GEdge*> > _temp;
+ std::set<myBundle<GEdge*> > _bundles;
+ std::map<MVertex*, GVertex*> _existingVertices;
+ std::map<GEdge*, std::set<GVertex*> > _topology;
+
+ // create an inverse dictionnary for existing edges
+ for(GModel::viter it = gm->firstVertex(); it != gm->lastVertex(); it++) {
+ _existingVertices[(*it)->mesh_vertices[0]] = *it;
+ }
+
+ // printf("%ud mesh vertices are already classified\n",_existingVertices.size());
+
+
+ for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++) {
+ for (int i=0;i<(*it)->lines.size();i++){
+ MLine *e = (*it)->lines[i];
+ for (int j=0;j<2;j++){
+ MVertex* v = e->getVertex(j);
+ GVertex *gv = _existingVertices[v];
+ if (gv)
+ _topology[*it].insert(gv);
+ else
+ _temp[v].insert(*it);
+ }
+ }
+ }
+
+ // printf("%d new mesh vertices \n",_temp.size());
+
+ // create unique instances
+ for (std::map<MVertex*, myBundle<GEdge*> >::iterator it = _temp.begin(); it != _temp.end() ; it++){
+ _bundles.insert (it->second);
+ }
+
+ // printf("%d bundles\n",_bundles.size());
+
+
+ // create discrete edges
+
+ std::map <myBundle<GEdge*> , GVertex *> _e2v;
+ {
+ std::set<myBundle<GEdge*> >::iterator it = _bundles.begin();
+ for (; it != _bundles.end(); ++it) {
+ if (it->stuff.size() > 1) {
+ discreteVertex *dv = new discreteVertex ( gm , NEWPOINT());
+ _e2v [*it] = dv;
+ for (std::set<GEdge*>::iterator it2 = it->stuff.begin(); it2 != it->stuff.end() ; ++it2){
+ _topology[*it2].insert(dv);
+ }
+ }
+ }
+ }
+
+ // create the 0D mesh
+ {
+ std::map<MVertex*, myBundle<GEdge*> > :: iterator it = _temp.begin();
+ for (; it != _temp.end(); ++it) {
+ if (it->second.stuff.size() > 1) {
+ MVertex *v = it->first;
+ GVertex *gv = _e2v [it->second];
+ gv->mesh_vertices.push_back(v);
+ }
+ }
+ }
+
+ // create Edge 2 Vertex topology
+ {
+ std::map<GEdge*, std::set<GVertex*> >::iterator it = _topology.begin();
+ for ( ; it != _topology.end() ; ++it){
+ std::list<GVertex*> l ; l.insert (l.begin(), it->second.begin(), it->second.end());
+ if (l.size() == 0){
+ }
+ else if (l.size() <= 2){
+ std::list<GVertex*>::iterator it2 = l.begin();
+ it->first->setBeginVertex(*it2);
+ if (l.size() == 2)++it2;
+ it->first->setEndVertex(*it2);
+ }
+ else {
+ Msg::Error ("FIXME : create simply connected edges in CreateTopology");
+ }
+
+ for (std::list<GVertex*>::iterator it2 = l.begin() ; it2 != l.end() ; ++it2)(*it2)->addEdge(it->first);
+ }
+ }
+
+
+ // NICE :-)
+}
+
+void createTopologyFromMesh2D ( GModel *gm ) {
+
+ std::map<MEdge, myBundle<GFace*>, Less_Edge > _temp;
+ std::set<myBundle<GFace*> > _bundles;
+ std::map<MEdge, GEdge*, Less_Edge> _existingEdges;
+ std::map<GFace*, std::set<GEdge*> > _topology;
+
+ // create an inverse dictionnary for existing edges
+ for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++) {
+ for (int i = 0; i < (*it)->lines.size(); i++)_existingEdges[(*it)->lines[i]->getEdge(0)] = *it;
+ }
+
+ // printf("%d mesh edges aere already classified\n",_existingEdges.size());
+
+ // create inverse dictionary for all other edges
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); it++) {
+ for (int i=0;i<(*it)->getNumMeshElements();i++){
+ MElement *e = (*it)->getMeshElement(i);
+ for (int j=0;j<e->getNumEdges();j++){
+ MEdge ed = e->getEdge(j);
+ GEdge *ge = _existingEdges[ed];
+ if (ge)
+ _topology[*it].insert(ge);
+ else
+ _temp[ed].insert(*it);
+ }
+ }
+ }
+
+ // create unique instances
+ for (std::map<MEdge, myBundle<GFace*>, Less_Face >::iterator it = _temp.begin(); it != _temp.end() ; it++){
+ _bundles.insert (it->second);
+ }
+
+ // create discrete edges
+ std::map <myBundle<GFace*> , GEdge *> _f2e;
+ {
+ std::set<myBundle<GFace*> >::iterator it = _bundles.begin();
+ for (; it != _bundles.end(); ++it) {
+ // it->print();
+ if (it->stuff.size() > 1){
+ printf("creation of a new discrete edge !\n");
+ discreteEdge *de = new discreteEdge ( gm , NEWREG(), NULL, NULL);
+ _f2e [*it] = de;
+ for (std::set<GFace*>::iterator it2 = it->stuff.begin(); it2 != it->stuff.end();++it2)
+ _topology[*it2].insert(de);
+ }
+ }
+ }
+
+ // create the 1D mesh
+ {
+ std::map<MEdge, myBundle<GFace*>, Less_Edge > :: iterator it = _temp.begin();
+ for (; it != _temp.end(); ++it) {
+ if (it->second.stuff.size() > 1){
+ MEdge e = it->first;
+ GEdge *ge = _f2e [it->second];
+ // not yet bounded
+ MLine *l = new MLine (e.getVertex(0),e.getVertex(1));
+ ge->lines.push_back(l);
+ for (int i=0;i<e.getNumVertices ();i++)e.getVertex(i)->setEntity(ge);
+ }
+ }
+ }
+
+ // create Face 2 Edge topology
+ {
+ std::map<GFace*, std::set<GEdge*> >::iterator it = _topology.begin();
+ for ( ; it != _topology.end() ; ++it){
+ std::list<GEdge*> l ; l.insert (l.begin(), it->second.begin(), it->second.end());
+ it->first->set(l);
+ // printf("Face %d has %d edges\n",it->first->tag(), l.size());
+ for (std::list<GEdge*>::iterator it2 = l.begin() ; it2 != l.end() ; ++it2)(*it2)->addFace(it->first);
+ }
+ }
+
+
+ // NICE :-)
+}
+
+/// ----------------------------------------------------------------
+/// --------------- 3D -----------------------------------
+/// ----------------------------------------------------------------
+
+void createTopologyFromMesh3D ( GModel *gm ) {
+
+ std::map<MFace, std::pair<GRegion*,GRegion*>, Less_Face > _temp;
+ std::set<std::pair<GRegion*,GRegion*> > _pairs;
+ std::map<MFace, GFace*, Less_Face> _existingFaces;
+ std::map<GRegion*, std::set<GFace*> > _topology;
+
+ // create an inverse dictionnary for existing faces
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); it++) {
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++)_existingFaces[(*it)->triangles[i]->getFace(0)] = *it;
+ for (unsigned int i = 0; i < (*it)->quadrangles.size(); i++)_existingFaces[(*it)->quadrangles[i]->getFace(0)] = *it;
+ }
+
+ // printf("%d mesh faces aere already classified\n",_existingFaces.size());
+
+ // create inverse dictionary for all other faces
+ for(GModel::riter it = gm->firstRegion(); it != gm->lastRegion(); it++) {
+ for (int i=0;i<(*it)->getNumMeshElements();i++){
+ MElement *e = (*it)->getMeshElement(i);
+ for (int j=0;j<e->getNumFaces();j++){
+ MFace f = e->getFace(j);
+ GFace *gf = _existingFaces[f];
+ if (gf){
+ _topology[*it].insert(gf);
+ }
+ else {
+ std::map<MFace, std::pair<GRegion*,GRegion*>, Less_Face >::iterator itf = _temp.find(f);
+ if (itf == _temp.end()){
+ _temp[f] = std::make_pair ( (GRegion*)NULL, *it) ;
+ }
+ else {
+ itf->second.first = *it ;
+ }
+ }
+ }
+ }
+ }
+
+ // printf("%d new mesh faces \n",_temp.size());
+
+ // create unique instances
+ for (std::map<MFace, std::pair<GRegion*,GRegion*>, Less_Face >::iterator it = _temp.begin(); it != _temp.end() ; it++){
+ _pairs.insert (std::make_pair (std::min (it->second.first, it->second.second),
+ std::max (it->second.first, it->second.second)));
+ }
+
+ // create discrete faces
+
+ // printf("%d pairs of regions exist to create new GFaces \n",_pairs.size());
+
+ std::map <std::pair<GRegion*,GRegion*> , GFace *> _r2f;
+ {
+ std::set<std::pair<GRegion*,GRegion*> >::iterator it = _pairs.begin();
+ for (; it != _pairs.end(); ++it) {
+ if (it->first != it->second) {
+ printf("creating a new discrete face between %p and %p\n",it->first, it->second);
+ discreteFace *df = new discreteFace ( gm , NEWREG());
+ _r2f [*it] = df;
+ if (it->first)_topology[it->first].insert(df);
+ if (it->second)_topology[it->second].insert(df);
+ }
+ }
+ }
+
+ // add mesh faces in newly created GFaces
+ {
+ std::map<MFace, std::pair<GRegion*,GRegion*>, Less_Face > :: iterator it = _temp.begin();
+ for (; it != _temp.end(); ++it) {
+ if (it->second.first != it->second.second){
+ MFace f = it->first;
+ GFace *gf = _r2f [it->second];
+ if (f.getNumVertices () == 3){
+ MTriangle *t = new MTriangle (f.getVertex(0),f.getVertex(1),f.getVertex(2));
+ gf->triangles.push_back(t);
+ }
+ else if (f.getNumVertices () == 4){
+ MQuadrangle *q = new MQuadrangle (f.getVertex(0),f.getVertex(1),f.getVertex(2),f.getVertex(3));
+ gf->quadrangles.push_back(q);
+ }
+ for (int i=0;i<f.getNumVertices ();i++)f.getVertex(i)->setEntity(gf);
+ }
+ }
+ }
+
+ // create Regions 2 Faces topology
+ {
+ std::map<GRegion*, std::set<GFace*> >::iterator it = _topology.begin();
+ for ( ; it != _topology.end() ; ++it){
+ std::list<GFace*> l ; l.insert (l.begin(), it->second.begin(), it->second.end());
+ it->first->set(l);
+ // printf("Region %d has %d adjacent faces\n",it->first->tag(), l.size());
+ for (std::list<GFace*>::iterator it2 = l.begin() ; it2 != l.end() ; ++it2)(*it2)->addRegion(it->first);
+ }
+ }
+
+ // NICE :-)
+
+}
+
+void createTopologyFromMeshNew ( GModel *gm ) {
+ const int dim = gm->getDim ();
+ if (dim >= 3) createTopologyFromMesh3D (gm);
+ if (dim >= 2) createTopologyFromMesh2D (gm);
+ if (dim >= 1) createTopologyFromMesh1D (gm);
+ Msg::Info("createTopologyFromMeshNew (%d regions %d faces and %d edges)",gm->getNumRegions(), gm->getNumFaces(), gm->getNumEdges());
+
+}
diff --git a/Geo/GModelCreateTopologyFromMesh.h b/Geo/GModelCreateTopologyFromMesh.h
new file mode 100644
index 0000000000000000000000000000000000000000..79a8ecf4f948a93ddff7f7394b552eda4492d664
--- /dev/null
+++ b/Geo/GModelCreateTopologyFromMesh.h
@@ -0,0 +1,5 @@
+#ifndef _CREATE_TOPOLOGY_FROM_MESH
+#define _CREATE_TOPOLOGY_FROM_MESH
+class GModel;
+void createTopologyFromMeshNew ( GModel *gm ) ;
+#endif
diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index d34eb6120103bf8b4e44837abbb722af1aec62f9..39f64b3ee882a278bdc6f15dc9343e3fbc197ee4 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -42,7 +42,6 @@ int GModel::readGEO(const std::string &name)
{
ParseFile(name, true);
int imported = GModel::current()->importGEOInternals();
-
return imported;
}
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index 04fe00a7ac6236e0a8fdf2f4a4964a8327094897..2e715cde930d28b465fba1f56a1641c0d3155d51 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/Geo/GModelIO_MSH.cpp
@@ -23,6 +23,7 @@
#include "discreteEdge.h"
#include "discreteFace.h"
#include "discreteRegion.h"
+#include "GModelCreateTopologyFromMesh.h"
static int readMSHPhysicals(FILE *fp, GEntity *ge)
{
@@ -179,6 +180,7 @@ int GModel::readMSH(const std::string &name)
double version = 0.;
bool binary = false, swap = false, postpro = false;
+ bool topology = false;
int minVertex = 0;
std::map<int, std::vector<MElement*> > elements[11];
@@ -234,6 +236,7 @@ int GModel::readMSH(const std::string &name)
// $Entities section
else if(!strncmp(&str[1], "Entities", 8)) {
readMSHEntities(fp, this);
+ topology = true;
}
// $Nodes section
@@ -444,6 +447,8 @@ int GModel::readMSH(const std::string &name)
} while(str[0] != '$');
}
+
+
// store the elements in their associated elementary entity. If the
// entity does not exist, create a new (discrete) one.
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
@@ -458,6 +463,9 @@ int GModel::readMSH(const std::string &name)
else
_storeVerticesInEntities(_vertexMapCache);
+ // if no topology is given, create one
+ if (!topology)createTopologyFromMeshNew (this);
+
_createGeometryOfDiscreteEntities() ;
for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
diff --git a/Geo/GModelIO_MSH2.cpp b/Geo/GModelIO_MSH2.cpp
index 2a6521750978b5e9ac697f75899d3f5022f52525..bcfae5b10618196feb1864e8bee3c8ec0b4531d5 100644
--- a/Geo/GModelIO_MSH2.cpp
+++ b/Geo/GModelIO_MSH2.cpp
@@ -26,6 +26,7 @@
#include "GmshMessage.h"
#include "Context.h"
#include "OS.h"
+#include "GModelCreateTopologyFromMesh.h"
#define FAST_ELEMENTS 1
@@ -687,6 +688,9 @@ int GModel::_readMSH2(const std::string &name)
else
_storeVerticesInEntities(vertexMap);
+ // if no topology is given, create one
+ createTopologyFromMeshNew (this);
+
// store the physical tags
for(int i = 0; i < 4; i++)
_storePhysicalTagsInEntities(i, physicals[i]);
diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index 6b1e41365fb0111faa5e564870a1e74e1de6750c..92295867ee1b51acc402d65539eb9771d4414900 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -741,6 +741,7 @@ GRegion* OCC_Internals::addRegionToModel(GModel *model, TopoDS_Solid region)
GRegion *gr = getOCCRegionByNativePtr(model, region);
if(gr) return gr;
+ // FIXME THE PREVIOUS IMPLEMENTATION WAS BETTER FOR SOME USERS :-)
buildShapeFromLists(region);
model->destroy();
buildLists();
diff --git a/Geo/GModelIO_STL.cpp b/Geo/GModelIO_STL.cpp
index a684d6e5dc6237285fe2f3a5e3bd20d3cc3bc02f..50abbcb315fe1bd3f6b681a7b37f5a717219a71f 100644
--- a/Geo/GModelIO_STL.cpp
+++ b/Geo/GModelIO_STL.cpp
@@ -12,6 +12,7 @@
#include "MVertexRTree.h"
#include "discreteFace.h"
#include "StringUtils.h"
+#include "GModelCreateTopologyFromMesh.h"
int GModel::readSTL(const std::string &name, double tolerance)
{
@@ -176,6 +177,11 @@ int GModel::readSTL(const std::string &name, double tolerance)
_associateEntityWithMeshVertices();
_storeVerticesInEntities(vertices); // will delete unused vertices
+ // if no topology is given, create one
+ createTopologyFromMeshNew (this);
+
+ _createGeometryOfDiscreteEntities() ;
+
fclose(fp);
return 1;
}
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 4f14d5c3d21043f1d4d323b12610726c2752f4ca..be03cb69193ff3fc55bba0e2126f055bf9163ac4 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -915,6 +915,7 @@ void RecombineMesh(GModel *m)
void GenerateMesh(GModel *m, int ask)
{
+
// ProfilerStart("gmsh.prof");
if(CTX::instance()->lock) {
Msg::Info("I'm busy! Ask me that later...");