From a5cd2b175af83f4e3a9b9f546a42a5c0098f87a5 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sun, 16 Nov 2008 08:39:01 +0000
Subject: [PATCH] j lechelle
---
Geo/GModelIO_Mesh.cpp | 105 ++++++++++++++++++++++--------------------
Geo/MElement.h | 52 +++++++++++++++++++--
2 files changed, 104 insertions(+), 53 deletions(-)
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 54db9b53be..93d5d3f083 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -2163,7 +2163,7 @@ int GModel::readDIFF(const std::string &name)
}
// FIXME: todo
-#if 0
+#if 0
char str[256] = "XXX";
std::map<int, std::vector<MElement*> > elements[8];
std::map<int, std::map<int, std::string> > physicals[4];
@@ -2221,7 +2221,7 @@ int GModel::readDIFF(const std::string &name)
Msg::Info("several_subdomains %x", several_subdomains);
int nbi;
- int* bi;
+ std::vector<int> bi;
if(!fgets(str, sizeof(str), fp) || feof(fp)) return 0;
while(strstr(str, "Boundary indicators:") == NULL){
if(!fgets(str, sizeof(str), fp) || feof(fp))
@@ -2229,8 +2229,8 @@ int GModel::readDIFF(const std::string &name)
}
if(sscanf(str, "%d %*s %*s", &nbi) != 1) return 0;
Msg::Info("nbi %d", nbi);
- if(nbi != 0)
- bi = new int[nbi];
+ if(nbi != 0)
+ bi.resize(nbi);
std::string format_read_bi = "%*d %*s %*s";
for(int i = 0; i < nbi; i++){
if(format_read_bi[format_read_bi.size()-1] == 'd') {
@@ -2239,7 +2239,7 @@ int GModel::readDIFF(const std::string &name)
}
else
format_read_bi += " %d";
- if(sscanf(str, format_read_bi.c_str(), bi + i) != 1) return 0;
+ if(sscanf(str, format_read_bi.c_str(), &bi[i]) != 1) return 0;
Msg::Info("bi[%d]=%d", i, bi[i]);
}
@@ -2251,16 +2251,17 @@ int GModel::readDIFF(const std::string &name)
vertexMap.clear();
int minVertex = numVertices + 1, maxVertex = -1;
int num;
- int **elementary;
- elementary = new int*[numVertices];
+ std::vector<std::vector<int> > elementary(numVertices);
Msg::ResetProgressMeter();
for(int i = 0; i < numVertices; i++){
- elementary[i] = new int[nbi + 1];
if(!fgets(str, sizeof(str), fp)) return 0;
double xyz[3];
+ int tmp;
if(sscanf(str, "%d ( %lf , %lf , %lf ) [%d]", &num, &xyz[0], &xyz[1], &xyz[2],
- &elementary[i][0]) != 5) return 0;
+ &tmp) != 5) return 0;
+ elementary[i].resize(tmp + 1);
+ elementary[i][0] = tmp;
minVertex = std::min(minVertex, num);
maxVertex = std::max(maxVertex, num);
if(vertexMap.count(num))
@@ -2303,49 +2304,55 @@ int GModel::readDIFF(const std::string &name)
break;
}
- int material[numElements];
- int ElementsNodes[numElements][numVerticesPerElement];
+ std::vector<int> material(numElements);
+ std::vector<std::vector<int> > ElementsNodes(numElements);
+ for(int i = 0; i < numVertices; i++){
+ ElementsNodes[i].resize(numVerticesPerElement);
+ }
+ char eleType[20]="";
Msg::ResetProgressMeter();
for(int i = 1; i <= numElements; i++){
- if(!fgets(str, sizeof(str), fp)) return 0;
- int num, type, physical = 0, partition = 0;
- int indices[60];
- if(numVerticesPerElement == 10){
- if(sscanf(str, "%d %*s %d %d %d %d %d %d %d %d %d %d %d", &num, &material[i - 1],
- &ElementsNodes[i - 1][1], &ElementsNodes[i - 1][0],
- &ElementsNodes[i - 1][2], &ElementsNodes[i - 1][3],
- &ElementsNodes[i - 1][4], &ElementsNodes[i - 1][6],
- &ElementsNodes[i - 1][5], &ElementsNodes[i - 1][9],
- &ElementsNodes[i - 1][7], &ElementsNodes[i - 1][8]) != 12) return 0;
- Msg::Info("%d %d %d %d %d %d %d %d %d %d %d %d", i, material[i - 1],
- ElementsNodes[i - 1][0], ElementsNodes[i - 1][1], ElementsNodes[i - 1][2],
- ElementsNodes[i - 1][3], ElementsNodes[i - 1][4], ElementsNodes[i - 1][5],
- ElementsNodes[i - 1][6], ElementsNodes[i - 1][7], ElementsNodes[i - 1][8],
- ElementsNodes[i - 1][9]);
- type = MSH_TET_10;
- }
- else {
- if(sscanf(str,"%d %*s %d %d %d %d %d", &num, &material[i - 1],
- &ElementsNodes[i - 1][1], &ElementsNodes[i - 1][0], &ElementsNodes[i - 1][2],
- &ElementsNodes[i - 1][3]) != 6) return 0;
- Msg::Info("%d %d %d %d %d %d", i, material[i - 1], ElementsNodes[i - 1][0],
- ElementsNodes[i - 1][1], ElementsNodes[i - 1][2], ElementsNodes[i - 1][3]);
- type = MSH_TET_4;
+ if(!fgets(str, sizeof(str), fp)) return 0;
+ int num, type, physical = 0, partition = 0;
+ int indices[60];
+ if(sscanf(str, "%*d %s %d", eleType, &material[i-1])!=2) return 0;
+ std::string format_read_vertices = "%*d %*s %*d";
+ for(int k = 0; k < numVerticesPerElement; k++){
+ if(format_read_vertices[format_read_vertices.size()-2] != '*') {
+ format_read_vertices[format_read_vertices.size()-1] = '*';
+ format_read_vertices += "d %d";
+ }
+ else
+ format_read_vertices += " %d";
+ int k2;
+ if(eleType=="ElmT10n3D"){
+ static int map[10]={1, 0, 2, 3, 4, 6, 5, 9, 7, 8};
+ k2=map[k];
+ type= MSH_TET_10;
}
- for(int j=0;j<numVerticesPerElement;j++)
- indices[j] = ElementsNodes[i - 1][j];
- std::vector<MVertex*> vertices;
- if(vertexVector.size()){
- if(!getVertices(numVerticesPerElement, indices, vertexVector, vertices)) return 0;
- }
- else{
- if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices)) return 0;
- }
- createElementMSH(this, 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, "Reading elements");
+ else if(eleType=="ElmT4n3D"){
+ static int map[4]={1, 0, 2, 3};
+ k2=map[k];
+ type= MSH_TET_4;
+ }
+ else
+ return 0;
+ if(sscanf(str, format_read_vertices.c_str(), &ElementsNodes[i-1][k2]) != 1) return 0;
+ }
+ for(int j=0;j<numVerticesPerElement;j++)
+ indices[j] = ElementsNodes[i - 1][j];
+ std::vector<MVertex*> vertices;
+ if(vertexVector.size()){
+ if(!getVertices(numVerticesPerElement, indices, vertexVector, vertices)) return 0;
+ }
+ else{
+ if(!getVertices(numVerticesPerElement, indices, vertexMap, vertices)) return 0;
+ }
+ createElementMSH(this, 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, "Reading elements");
}
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 0935f163a3..adc49b1b0e 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -92,9 +92,8 @@ class MElement
// get the vertex using MED ordering
virtual MVertex *getVertexMED(int num){ return getVertex(num); }
- // get the vertex using DIFF ordering (at least for tetrahedra it's
- // the same as in the MED format)
- virtual MVertex *getVertexDIFF(int num){ return getVertexMED(num); }
+ // get the vertex using DIFF ordering
+ virtual MVertex *getVertexDIFF(int num){ return getVertex(num); }
// get the number of vertices associated with edges, faces and
// volumes (nonzero only for higher order elements)
@@ -572,6 +571,7 @@ class MTriangle : public MElement {
virtual int getTypeForVTK() const { return 5; }
virtual const char *getStringForPOS() const { return "ST"; }
virtual const char *getStringForBDF() const { return "CTRIA3"; }
+ virtual const char *getStringForDIFF() const { return "ElmT3n2D"; }
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
@@ -679,6 +679,7 @@ class MTriangle6 : public MTriangle {
//virtual int getTypeForVTK() const { return 22; }
virtual const char *getStringForPOS() const { return "ST2"; }
virtual const char *getStringForBDF() const { return "CTRIA6"; }
+ virtual const char *getStringForDIFF() const { return "ElmT6n2D"; }
virtual void revert()
{
MVertex *tmp;
@@ -832,6 +833,11 @@ class MQuadrangle : public MElement {
static const int map[4] = {0, 3, 2, 1};
return getVertex(map[num]);
}
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ static const int map[4] = {0, 1, 3, 2};
+ return getVertex(map[num]);
+ }
virtual int getNumEdges(){ return 4; }
virtual MEdge getEdge(int num)
{
@@ -869,6 +875,7 @@ class MQuadrangle : public MElement {
virtual int getTypeForVTK() const { return 9; }
virtual const char *getStringForPOS() const { return "SQ"; }
virtual const char *getStringForBDF() const { return "CQUAD4"; }
+ virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
virtual void revert()
{
MVertex *tmp = _v[1]; _v[1] = _v[3]; _v[3] = tmp;
@@ -950,6 +957,11 @@ class MQuadrangle8 : public MQuadrangle {
static const int map[8] = {0, 3, 2, 1, 7, 6, 5, 4};
return getVertex(map[num]);
}
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ static const int map[8] = {0, 1, 3, 2, 4, 7, 5, 6};
+ return getVertex(map[num]);
+ }
virtual int getNumEdgeVertices() const { return 4; }
virtual int getNumEdgesRep(){ return 8; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -990,6 +1002,7 @@ class MQuadrangle8 : public MQuadrangle {
virtual int getTypeForUNV() const { return 95; } // shell parabolic quadrilateral
//virtual int getTypeForVTK() const { return 23; }
virtual const char *getStringForBDF() const { return "CQUAD8"; }
+ virtual const char *getStringForDIFF() const { return "ElmB8n2D"; }
virtual void revert()
{
MVertex *tmp;
@@ -1032,8 +1045,13 @@ class MQuadrangle9 : public MQuadrangle {
virtual int getPolynomialOrder() const { return 2; }
virtual int getNumVertices() const { return 9; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ static const int map[9] = {0, 2, 8, 6, 1, 5, 7, 3, 4};
+ return getVertex(map[num]);
+ }
virtual int getNumEdgeVertices() const { return 4; }
- virtual int getNumFaceVertices() const { return 1; }
+ virtual int getNumFaceVertices() const { return 6; }
virtual int getNumEdgesRep(){ return 8; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
@@ -1073,6 +1091,7 @@ class MQuadrangle9 : public MQuadrangle {
}
virtual int getTypeForMSH() const { return MSH_QUA_9; }
virtual const char *getStringForPOS() const { return "SQ2"; }
+ virtual const char *getStringForDIFF() const { return "ElmB9n2D"; }
virtual void revert()
{
MVertex *tmp;
@@ -1290,6 +1309,11 @@ class MTetrahedron10 : public MTetrahedron {
static const int map[10] = {0, 2, 1, 3, 6, 5, 4, 7, 8, 9};
return getVertex(map[num]);
}
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ static const int map[10] = {0, 1, 2, 3, 4, 5, 6, 7, 9, 8};
+ return getVertex(map[num]);
+ }
virtual int getNumEdgeVertices() const { return 6; }
virtual int getNumEdgesRep(){ return 12; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1584,6 +1608,11 @@ class MHexahedron : public MElement {
static const int map[8] = {0, 3, 2, 1, 4, 7, 6, 5};
return getVertex(map[num]);
}
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ static const int map[8] = {2, 3, 7, 6, 0, 1, 5, 4};
+ return getVertex(map[num]);
+ }
virtual int getNumEdges(){ return 12; }
virtual MEdge getEdge(int num)
{
@@ -1633,6 +1662,7 @@ class MHexahedron : public MElement {
virtual int getTypeForVTK() const { return 12; }
virtual const char *getStringForPOS() const { return "SH"; }
virtual const char *getStringForBDF() const { return "CHEXA"; }
+ virtual const char *getStringForDIFF() const { return "ElmB8n3D"; }
virtual void revert()
{
MVertex *tmp;
@@ -1779,6 +1809,12 @@ class MHexahedron20 : public MHexahedron {
8, 17, 19, 18, 16, 10, 15, 14, 12};
return getVertex(map[num]);
}
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ static const int map[20] = {2, 3, 7, 6, 0, 1, 5, 4, 9, 18, 12,
+ 19, 14, 11, 15, 13, 8, 16, 17, 10};
+ return getVertex(map[num]);
+ }
virtual int getNumEdgeVertices() const { return 12; }
virtual int getNumEdgesRep(){ return 24; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
@@ -1841,6 +1877,7 @@ class MHexahedron20 : public MHexahedron {
virtual int getTypeForUNV() const { return 116; } // solid parabolic brick
//virtual int getTypeForVTK() const { return 25; }
virtual const char *getStringForBDF() const { return "CHEXA"; }
+ virtual const char *getStringForDIFF() const { return "ElmB20n3D"; }
virtual void revert()
{
MVertex *tmp;
@@ -1900,6 +1937,12 @@ class MHexahedron27 : public MHexahedron {
virtual int getPolynomialOrder() const { return 2; }
virtual int getNumVertices() const { return 27; }
virtual MVertex *getVertex(int num){ return num < 8 ? _v[num] : _vs[num - 8]; }
+ virtual MVertex *getVertexDIFF(int num)
+ {
+ static const int map[27] = {6, 8, 26, 24, 0, 2, 20, 18, 7, 15, 3, 17, 5, 25,
+ 23, 21, 1, 9, 11, 19, 16, 4, 12, 14, 22, 10, 13};
+ return getVertex(map[num]);
+ }
virtual int getNumEdgeVertices() const { return 12; }
virtual int getNumFaceVertices() const { return 6; }
virtual int getNumVolumeVertices() const { return 1; }
@@ -1969,6 +2012,7 @@ class MHexahedron27 : public MHexahedron {
}
virtual int getTypeForMSH() const { return MSH_HEX_27; }
virtual const char *getStringForPOS() const { return "SH2"; }
+ virtual const char *getStringForDIFF() const { return "ElmB27n3D"; }
virtual void revert()
{
MVertex *tmp;
--
GitLab