diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 8d2ece39a7766075d83c62334ae25e25e4d1bf02..e53dbd37a4e61955c08d7ae4adfd2dce1b640e26 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1237,10 +1237,10 @@ static bool getVertices(int num, int *indices, std::vector<MVertex*> &vec,
return true;
}
-MElement *GModel::createElementMSH(int num, int typeMSH, int physical,
- int reg, int part, std::vector<MVertex*> &v,
- std::map<int, std::vector<MElement*> > elements[10],
- std::map<int, std::map<int, std::string> > physicals[4])
+MElement *GModel::_createElementMSH(int num, int typeMSH, int physical,
+ int reg, int part, std::vector<MVertex*> &v,
+ std::map<int, std::vector<MElement*> > elements[10],
+ std::map<int, std::map<int, std::string> > physicals[4])
{
if(CTX::instance()->mesh.switchElementTags) {
int tmp = reg;
@@ -1355,8 +1355,8 @@ GModel *GModel::createGModel(std::map<int, MVertex*> &vertexMap,
}
}
- gm->createElementMSH(num, elementType[i], physical[i], elementary[i],
- partition[i], vertices, elements, physicals);
+ gm->_createElementMSH(num, elementType[i], physical[i], elementary[i],
+ partition[i], vertices, elements, physicals);
}
// store the elements in their associated elementary entity. If the
@@ -2361,6 +2361,7 @@ GEdge* GModel::addCompoundEdge(std::vector<GEdge*> edges, int num){
return gec;
}
+
GFace* GModel::addCompoundFace(std::vector<GFace*> faces, int param, int split, int num)
{
#if defined(HAVE_SOLVER)
@@ -2500,7 +2501,6 @@ GFace* GModel::addFace (std::vector<GEdge *> edges,
return 0;
}
-
GFace* GModel::addPlanarFace (std::vector<std::vector<GEdge *> > edges){
if(_factory)
return _factory->addPlanarFace(this, edges);
@@ -2602,8 +2602,9 @@ static void computeDuplicates(GModel *model,
std::map<GVertex*, GVertex*> &Duplicates2Unique,
const double &eps)
{
- // FIXME: currently we use a greedy algorithm in n^2 (using a
- // kd-tree: cf. MVertexPositionSet)
+ // FIXME: currently we use a greedy algorithm in n^2 (using a kd-tree:
+ // cf. MVertexPositionSet)
+
// FIXME: add option to remove orphaned entities after duplicate check
std::list<GVertex*> v;
v.insert(v.begin(), model->firstVertex(), model->lastVertex());
diff --git a/Geo/GModel.h b/Geo/GModel.h
index aee48a484631483e11c2f7130a35528a9d31b29e..ea9b32a5e904fc0641b933b2cdf41b763d69b39c 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -123,6 +123,13 @@ class GModel
void _storePhysicalTagsInEntities(int dim,
std::map<int, std::map<int, std::string> > &map);
+ // utility function to create a mesh element in the I/O routines (don't use in
+ // new code; it will be removed eventually)
+ MElement *_createElementMSH(int num, int typeMSH, int physical,
+ int reg, int part, std::vector<MVertex*> &v,
+ std::map<int, std::vector<MElement*> > elements[10],
+ std::map<int, std::map<int, std::string> > physicals[4]);
+
// entity that is currently being meshed (used for error reporting)
GEntity *_currentMeshEntity;
@@ -488,15 +495,9 @@ class GModel
// if cutElem is set to false, split the model without cutting the elements
GModel *buildCutGModel(gLevelset *ls, bool cutElem=true, bool saveTri=false);
- // utility function to create a mesh element in the I/O routines
- MElement *createElementMSH(int num, int typeMSH, int physical,
- int reg, int part, std::vector<MVertex*> &v,
- std::map<int, std::vector<MElement*> > elements[10],
- std::map<int, std::map<int, std::string> > physicals[4]);
-
- // 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
- // to the number of elements)
+ // 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 to the number
+ // of elements)
static GModel *createGModel(std::map<int, MVertex*> &vertexMap,
std::vector<int> &numElement,
std::vector<std::vector<int> > &vertexIndices,
@@ -505,9 +506,9 @@ 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
+ // 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);
diff --git a/Geo/GModelIO_DIFF.cpp b/Geo/GModelIO_DIFF.cpp
index 079fe25b159eb9361f50840fe70d80874ca45f3f..0e5310f061a506a835233e0b5e31642d38e42aad 100644
--- a/Geo/GModelIO_DIFF.cpp
+++ b/Geo/GModelIO_DIFF.cpp
@@ -294,8 +294,8 @@ int GModel::readDIFF(const std::string &name)
if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices))
return 0;
}
- createElementMSH(num, type, physical, elementary[i-1][1], partition,
- vertices, elements, physicals);
+ _createElementMSH(num, type, physical, elementary[i-1][1], partition,
+ vertices, elements, physicals);
// trouble if elementary[i-1][0]>1 nodal post-processing needed ?
if(numElements > 100000)
Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
diff --git a/Geo/GModelIO_MSH.cpp b/Geo/GModelIO_MSH.cpp
index 0401e647b52f451c9eb6f24ae7b475b587620dfe..3e949fd2fb63b20c1bc0a7acac77a77777dccb85 100644
--- a/Geo/GModelIO_MSH.cpp
+++ b/Geo/GModelIO_MSH.cpp
@@ -14,14 +14,345 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MPyramid.h"
+#include "StringUtils.h"
+
+static bool getVertices(int num, const std::vector<int> &indices,
+ std::map<int, MVertex*> &map,
+ std::vector<MVertex*> &vertices)
+{
+ for(int i = 0; i < num; i++){
+ if(!map.count(indices[i])){
+ Msg::Error("Wrong vertex index %d", indices[i]);
+ return false;
+ }
+ else
+ vertices.push_back(map[indices[i]]);
+ }
+ return true;
+}
+
+static bool getVertices(int num, const std::vector<int> &indices,
+ std::vector<MVertex*> &vec,
+ std::vector<MVertex*> &vertices, int minVertex = 0)
+{
+ for(int i = 0; i < num; i++){
+ if(indices[i] < minVertex || indices[i] > (int)(vec.size() - 1 + minVertex)){
+ Msg::Error("Wrong vertex index %d", indices[i]);
+ return false;
+ }
+ else
+ vertices.push_back(vec[indices[i]]);
+ }
+ return true;
+}
int GModel::readMSH(const std::string &name)
{
+ FILE *fp = fopen(name.c_str(), "rb");
+ if(!fp){
+ Msg::Error("Unable to open file '%s'", name.c_str());
+ return 0;
+ }
+
+ char str[256] = "";
double version = 2.2;
- if(version < 3)
- return _readMSH2(name);
+ bool binary = false, swap = false, postpro = false;
+ int minVertex = 0;
+ std::map<int, std::vector<int> > entities[4];
+ std::map<int, std::vector<MElement*> > elements[10];
+
+ while(1) {
+
+ while(str[0] != '$'){
+ if(!fgets(str, sizeof(str), fp) || feof(fp))
+ break;
+ }
+
+ if(feof(fp))
+ break;
+
+ // $MeshFormat section
+ if(!strncmp(&str[1], "MeshFormat", 10)) {
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ int format, size;
+ if(sscanf(str, "%lf %d %d", &version, &format, &size) != 3) return 0;
+ if(version < 3.) return _readMSH2(name);
+ if(format){
+ binary = true;
+ Msg::Info("Mesh is in binary format");
+ int one;
+ if(fread(&one, sizeof(int), 1, fp) != 1) return 0;
+ if(one != 1){
+ swap = true;
+ Msg::Info("Swapping bytes from binary file");
+ }
+ }
+ }
+
+ // $PhysicalNames section
+ else if(!strncmp(&str[1], "PhysicalNames", 13)) {
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ int numNames;
+ if(sscanf(str, "%d", &numNames) != 1) return 0;
+ for(int i = 0; i < numNames; i++) {
+ int dim, num;
+ if(fscanf(fp, "%d", &dim) != 1) return 0;
+ if(fscanf(fp, "%d", &num) != 1) return 0;
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ std::string name = ExtractDoubleQuotedString(str, 256);
+ if(name.size()) setPhysicalName(name, dim, num);
+ }
+ }
+
+ // $Entities section
+ else if(!strncmp(&str[1], "Entities", 8)) {
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ int numEntities;
+ if(sscanf(str, "%d", &numEntities) != 1) return 0;
+ for(int i = 0; i < numEntities; i++) {
+ int num, dim, numPhysicals;
+ if(fscanf(fp, "%d %d %d", &num, &dim, &numPhysicals) != 3) return 0;
+ if(numPhysicals > 0){
+ std::vector<int> physicals(numPhysicals);
+ for(int j = 0; j < numPhysicals; j++){
+ if(fscanf(fp, "%d", &physicals[j]) != 1) return 0;
+ }
+ entities[dim][num] = physicals;
+ }
+ }
+ }
+
+ // $Nodes section
+ else if(!strncmp(&str[1], "Nodes", 5)) {
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ int numVertices;
+ if(sscanf(str, "%d", &numVertices) != 1) return 0;
+ Msg::Info("%d vertices", numVertices);
+ Msg::ResetProgressMeter();
+ std::map<int, MVertex*> vertexMap;
+ std::vector<MVertex*> vertexVector;
+ int maxVertex = -1;
+ minVertex = numVertices + 1;
+ for(int i = 0; i < numVertices; i++) {
+ int num, entity, dim;
+ double xyz[3];
+ MVertex *vertex = 0;
+ if(!binary){
+ if(fscanf(fp, "%d %lf %lf %lf %d %d", &num, &xyz[0], &xyz[1], &xyz[2],
+ &entity, &dim) != 6)
+ return 0;
+ }
+ else{
+ if(fread(&num, sizeof(int), 1, fp) != 1) return 0;
+ if(swap) SwapBytes((char*)&num, sizeof(int), 1);
+ if(fread(xyz, sizeof(double), 3, fp) != 3) return 0;
+ if(swap) SwapBytes((char*)xyz, sizeof(double), 3);
+ if(fread(&entity, sizeof(int), 1, fp) != 1) return 0;
+ if(swap) SwapBytes((char*)&entity, sizeof(int), 1);
+ if(fread(&dim, sizeof(int), 1, fp) != 1) return 0;
+ if(swap) SwapBytes((char*)&dim, sizeof(int), 1);
+ }
+ if(entity){
+ switch(dim){
+ case 0:
+ {
+ GVertex *gv = getVertexByTag(entity);
+ if (gv) gv->deleteMesh();
+ vertex = new MVertex(xyz[0], xyz[1], xyz[2], gv, num);
+ }
+ break;
+ case 1:
+ {
+ GEdge *ge = getEdgeByTag(entity);
+ double u;
+ if(!binary){
+ if (fscanf(fp, "%lf", &u) != 1) return 0;
+ }
+ else{
+ if(fread(&u, sizeof(double), 1, fp) != 1) return 0;
+ if(swap) SwapBytes((char*)&u, sizeof(double), 1);
+ }
+ vertex = new MEdgeVertex(xyz[0], xyz[1], xyz[2], ge, u, -1.0, num);
+ }
+ break;
+ case 2:
+ {
+ GFace *gf = getFaceByTag(entity);
+ double uv[2];
+ if(!binary){
+ if (fscanf(fp, "%lf %lf", &uv[0], &uv[1]) != 2) return 0;
+ }
+ else{
+ if(fread(uv, sizeof(double), 2, fp) != 2) return 0;
+ if(swap) SwapBytes((char*)uv, sizeof(double), 2);
+ }
+ vertex = new MFaceVertex(xyz[0], xyz[1], xyz[2], gf, uv[0], uv[1], num);
+ }
+ break;
+ case 3:
+ {
+ GRegion *gr = getRegionByTag(entity);
+ double uvw[3];
+ if(!binary){
+ if(fscanf(fp, "%lf %lf %lf", &uvw[0], &uvw[1], &uvw[2]) != 2) return 0;
+ }
+ else{
+ if(fread(uvw, sizeof(double), 3, fp) != 3) return 0;
+ if(swap) SwapBytes((char*)uvw, sizeof(double), 3);
+ }
+ vertex = new MVertex(xyz[0], xyz[1], xyz[2], gr, num);
+ }
+ break;
+ default:
+ Msg::Error("Wrong entity dimension for vertex %d", num);
+ return 0;
+ }
+ }
+ minVertex = std::min(minVertex, num);
+ maxVertex = std::max(maxVertex, num);
+ if(vertexMap.count(num))
+ Msg::Warning("Skipping duplicate vertex %d", num);
+ vertexMap[num] = vertex;
+ if(numVertices > 100000)
+ Msg::ProgressMeter(i + 1, numVertices, true, "Reading nodes");
+ }
+ // 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) ||
+ (minVertex == 0 && maxVertex == numVertices - 1))){
+ Msg::Info("Vertex numbering is dense");
+ vertexVector.resize(vertexMap.size() + 1);
+ if(minVertex == 1)
+ vertexVector[0] = 0;
+ else
+ vertexVector[numVertices] = 0;
+ std::map<int, MVertex*>::const_iterator it = vertexMap.begin();
+ for(; it != vertexMap.end(); ++it)
+ vertexVector[it->first] = it->second;
+ vertexMap.clear();
+ }
+
+ // cache the vertex indexing data
+ _vertexVectorCache = vertexVector;
+ _vertexMapCache = vertexMap;
+ }
+
+ // $Elements section
+ else if(!strncmp(&str[1], "Elements", 8)) {
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ int numElements;
+ if(sscanf(str, "%d", &numElements) != 1) return 0;
+ Msg::Info("%d elements", numElements);
+ Msg::ResetProgressMeter();
+ std::map<int, MElement*> elementMap;
+ for(int i = 0; i < numElements; i++) {
+ int num, type, numTags, numVertices;
+ if(!binary){
+ if(fscanf(fp, "%d %d %d", &num, &type, &numTags) != 3) return 0;
+ }
+ else{
+ if(fread(&num, sizeof(int), 1, fp) != 1) return 0;
+ if(swap) SwapBytes((char*)&num, sizeof(int), 1);
+ if(fread(&type, sizeof(int), 1, fp) != 1) return 0;
+ if(swap) SwapBytes((char*)&type, sizeof(int), 1);
+ if(fread(&numTags, sizeof(int), 1, fp) != 1) return 0;
+ if(swap) SwapBytes((char*)&numTags, sizeof(int), 1);
+ }
+ std::vector<int> tags;
+ if(numTags > 0){
+ tags.resize(numTags);
+ if(!binary){
+ for(int j = 0; j < numTags; j++){
+ if(fscanf(fp, "%d", &tags[j]) != 1) return 0;
+ }
+ }
+ else{
+ if(fread(&tags[0], sizeof(int), numTags, fp) != numTags) return 0;
+ if(swap) SwapBytes((char*)&tags[0], sizeof(int), numTags);
+ }
+ }
+ if(!(numVertices = MElement::getInfoMSH(type))) {
+ return 0;
+ }
+ std::vector<int> indices(numVertices);
+ if(!binary){
+ for(int j = 0; j < numVertices; j++)
+ if(fscanf(fp, "%d", &indices[j]) != 1) return 0;
+ }
+ else{
+ if(fread(&indices[0], sizeof(int), numVertices, fp) != numVertices) return 0;
+ if(swap) SwapBytes((char*)&indices[0], sizeof(int), numVertices);
+ }
+ std::vector<MVertex*> vertices;
+ if(_vertexVectorCache.size()){
+ if(!getVertices(numVertices, indices, _vertexVectorCache, vertices, minVertex))
+ return 0;
+ }
+ else{
+ if(!getVertices(numVertices, indices, _vertexMapCache, vertices))
+ return 0;
+ }
+ std::vector<short> ghosts;
+ MElementFactory factory;
+ MElement *element = factory.create(num, type, tags, vertices, elementMap,
+ ghosts);
+ elementMap[num] = element;
+ int part = element->getPartition();
+ if(part) getMeshPartitions().insert(part);
+ int elementary = tags.size() ? tags[0] : 0;
+ switch(element->getType()){
+ case TYPE_PNT: elements[0][elementary].push_back(element); break;
+ case TYPE_LIN: elements[1][elementary].push_back(element); break;
+ case TYPE_TRI: elements[2][elementary].push_back(element); break;
+ case TYPE_QUA: elements[3][elementary].push_back(element); break;
+ case TYPE_TET: elements[4][elementary].push_back(element); break;
+ case TYPE_HEX: elements[5][elementary].push_back(element); break;
+ case TYPE_PRI: elements[6][elementary].push_back(element); break;
+ case TYPE_PYR: elements[7][elementary].push_back(element); break;
+ }
+ for(unsigned int j = 0; j < ghosts.size(); j++)
+ _ghostCells.insert(std::pair<MElement*, short>(element, ghosts[j]));
+ if(numElements > 100000)
+ Msg::ProgressMeter(i + 1, numElements, true, "Reading elements");
+ }
+ // cache the element map
+ _elementMapCache = elementMap;
+ }
+
+ // Post-processing sections
+ else if(!strncmp(&str[1], "NodeData", 8) ||
+ !strncmp(&str[1], "ElementData", 11) ||
+ !strncmp(&str[1], "ElementNodeData", 15)) {
+ postpro = true;
+ break;
+ }
+
+ do {
+ if(!fgets(str, sizeof(str), fp) || feof(fp))
+ break;
+ } 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++)
+ _storeElementsInEntities(elements[i]);
+
+ // associate the correct geometrical entity with each mesh vertex
+ _associateEntityWithMeshVertices();
+
+ // store the vertices in their associated geometrical entity
+ if(_vertexVectorCache.size())
+ _storeVerticesInEntities(_vertexVectorCache);
+ else
+ _storeVerticesInEntities(_vertexMapCache);
+
+ // FIXME: store physicals in entities
+
+ fclose(fp);
- return 0;
+ return postpro ? 2 : 1;
}
static int getNumElementsMSH(GEntity *ge, bool saveAll, int saveSinglePartition)
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 87bc12f812bca3e194642ac3e00f0d1495aa7c34..02a359fb22c7f0b2ae0cdeae6d7af0eb59f173bc 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1398,6 +1398,32 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
}
}
+MElement *MElementFactory::create(int num, int type, const std::vector<int> &tags,
+ std::vector<MVertex*> &v,
+ std::map<int, MElement*> &elementCache,
+ std::vector<short> &ghosts)
+{
+ MElement *parent = 0;
+ int part = 0;
+ if(tags.size() > 2 && (type == MSH_PNT_SUB || type == MSH_LIN_SUB ||
+ type == MSH_TRI_SUB || type == MSH_TET_SUB)){
+ parent = elementCache[tags[1]];
+ if(tags[2]){ // num partitions
+ part = tags[3];
+ for(int i = 0; i < tags[2] - 1; i++)
+ ghosts.push_back(tags[4 + i]);
+ }
+ }
+ else if(tags.size() > 1){
+ if(tags[1]){ // num partitions
+ part = tags[2];
+ for(int i = 0; i < tags[1] - 1; i++)
+ ghosts.push_back(tags[3 + i]);
+ }
+ }
+ create(type, v, num, part, false, parent);
+}
+
void MElement::xyzTouvw(fullMatrix<double> *xu)
{
double _xyz[3] = {(*xu)(0,0),(*xu)(0,1),(*xu)(0,2)}, _uvw[3];
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 3ce5d69ce73b71c35a856fa1dbdf0c65d77c5b95..d1d33ad50cecffefe80033c02250379275b0c1e1 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -361,6 +361,10 @@ class MElementFactory{
public:
MElement *create(int type, std::vector<MVertex*> &v, int num=0, int part=0,
bool owner=false, MElement *parent=0, MElement *d1=0, MElement *d2=0);
+ MElement *create(int num, int type, const std::vector<int> &tags,
+ std::vector<MVertex*> &v,
+ std::map<int, MElement*> &elementCache,
+ std::vector<short> &ghosts);
};
// Traits of various elements based on the dimension. These generally define