From 8ad1324780220e4dc84c161006c0427453090124 Mon Sep 17 00:00:00 2001
From: Matti Pellika <matti.pellikka@tut.fi>
Date: Thu, 26 Apr 2012 09:13:25 +0000
Subject: [PATCH] New createGModel method.

---
 Geo/CellComplex.cpp   |  1 +
 Geo/GModel.cpp        |  8 ++++----
 Geo/GModel.h          |  7 +++++++
 Geo/GModelIO_Mesh.cpp | 46 ++++++++++++++++++++++++++++++++++++++-----
 4 files changed, 53 insertions(+), 9 deletions(-)

diff --git a/Geo/CellComplex.cpp b/Geo/CellComplex.cpp
index 1f34fe35d0..9359811aa0 100644
--- a/Geo/CellComplex.cpp
+++ b/Geo/CellComplex.cpp
@@ -433,6 +433,7 @@ int CellComplex::coreduceComplex(bool docombine, bool omit)
     for(unsigned int i = 0; i < newCells.size(); i++){
       insertCell(newCells.at(i));
     }
+
   }
 
   Msg::Debug(" %d volumes, %d faces, %d edges and %d vertices",
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 658660c04a..6bdde900fa 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -579,8 +579,8 @@ int GModel::adaptMesh(std::vector<int> technique, std::vector<simpleFunction<dou
   }
   //adapt only upper most dimension
   else{
-    
-    while(1) {  
+
+    while(1) {
       Msg::Info("-- adaptMesh ITER =%d ", ITER);
       std::vector<MElement*> elements;
 
@@ -2132,7 +2132,7 @@ GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split,
   if (param == 7) typ =  GFaceCompound::CONFORMAL_FE;
 
   GFaceCompound *gfc = new GFaceCompound(this, num, faces_comp, U0, typ, split);
-  
+
   //create old geo format for the compound face
   //necessary for boundary layers
   if(FindSurface(num)){
@@ -2171,7 +2171,7 @@ GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split,
      	}
        }
      }
-    
+
     Tree_Add(_geo_internals->Surfaces, &s);
   }
 
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 0c9b84f37d..7a890dc271 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -489,6 +489,13 @@ class GModel
                               std::vector<int> &elementary,
                               std::vector<int> &partition);
 
+  // create a GModel from newly created mesh elements
+  // (with their own newly created mesh vertices),
+  // and let element entities have given physical group tags
+  static GModel *createGModel
+    (std::map<int, std::vector<MElement*> > &entityToElementsMap,
+     std::map<int, std::vector<int> > &entityToPhysicalsMap);
+
   // for elements cut having new vertices
   void store(std::vector<MVertex*> &vertices, int dim,
             std::map<int, std::vector<MElement*> > &entityMap,
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 127492ebb7..741c601506 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -530,20 +530,20 @@ int GModel::readMSH(const std::string &name)
     MElement *e = createElementMSH(this, num, type, physical, elementary,
                                    partition, vertices, elements, physicals,
                                    own, p, doms[0], doms[1]);
-    
+
 #if (FAST_ELEMENTS==1)
 	  elems[num] = e;
 	  elemreg[num] = elementary;
 	  elemphy[num] = physical;
 #endif
-    
+
     for(unsigned int j = 0; j < ghosts.size(); j++)
       _ghostCells.insert(std::pair<MElement*, short>(e, ghosts[j]));
     if(numElements > 100000)
       Msg::ProgressMeter(i + 1, numElements, "Reading elements");
     delete [] indices;
         }
-        
+
 #if (FAST_ELEMENTS==1)
 	for(int i = 0; i < 10; i++)
 	  elements[i].clear();
@@ -3782,15 +3782,15 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
     nbVertices = (int)vertexIndices[i].size();
     indices = &vertexIndices[i][0];
     if(vertexVector.size()){
-      Msg::Error("Vertex not found aborting");
       if(!getVertices(nbVertices, indices, vertexVector, vertices)){
+        Msg::Error("Vertex not found aborting");
         delete gm;
         return 0;
       }
     }
     else{
-      Msg::Error("Vertex not found aborting");
       if(!getVertices(nbVertices, indices, vertexMap, vertices)){
+        Msg::Error("Vertex not found aborting");
         delete gm;
         return 0;
       }
@@ -3820,3 +3820,39 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
 
   return gm;
 }
+
+GModel *GModel::createGModel
+(std::map<int, std::vector<MElement*> > &entityToElementsMap,
+ std::map<int, std::vector<int> > &entityToPhysicalsMap)
+{
+  GModel* gm = new GModel();
+
+  std::map<int, MVertex*> vertexMap;
+  std::map<int, std::map<int, std::string> > physicals[4];
+  for(std::map<int, std::vector<MElement*> >::iterator it =
+        entityToElementsMap.begin(); it != entityToElementsMap.end();
+      it++) {
+    int entity = it->first;
+    for(unsigned int iE = 0; iE < it->second.size(); iE++) {
+      MElement* me = it->second[iE];
+      for(int iV = 0; iV < me->getNumVertices(); iV++) {
+        vertexMap[me->getVertex(iV)->getNum()] = me->getVertex(iV);
+      }
+      if(me->getPartition()) {
+        gm->getMeshPartitions().insert(me->getPartition());
+      }
+      std::vector<int> entityPhysicals = entityToPhysicalsMap[entity];
+      for(unsigned int i = 0; i < entityPhysicals.size(); i++) {
+        physicals[me->getDim()][entity][entityPhysicals[i]] = "unnamed";
+      }
+    }
+  }
+
+  gm->_storeElementsInEntities(entityToElementsMap);
+  gm->_associateEntityWithMeshVertices();
+  gm->_storeVerticesInEntities(vertexMap);
+  for(int i = 0; i < 4; i++)
+    gm->_storePhysicalTagsInEntities(i, physicals[i]);
+
+  return gm;
+}
-- 
GitLab