Skip to content
Snippets Groups Projects
Commit c3be78b0 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files 

- Use correct vertex ordering for second order simplices in UNV format
(thanks to Romuald Conty <conty@phimeca.com>)

- MLine2 -> MLine3
parent ce4e915b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/GModelIO.cpp
+ 71
− 18
View file @ c3be78b0
// $Id: GModelIO.cpp,v 1.43 2006-09-06 18:42:24 geuzaine Exp $ // $Id: GModelIO.cpp,v 1.44 2006-09-07 19:45:15 geuzaine Exp $
// //
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
// //
...@@ -236,7 +236,7 @@ static void createElementMSH(GModel *m, int num, int type, int physical, ...@@ -236,7 +236,7 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
break; break;
case LGN2: case LGN2:
elements[0][elementary].push_back elements[0][elementary].push_back
(new MLine2(v[n[0]], v[n[1]], v[n[2]], num, partition)); (new MLine3(v[n[0]], v[n[1]], v[n[2]], num, partition));
dim = 1; dim = 1;
break; break;
case TRI1: case TRI1:
...@@ -1186,24 +1186,77 @@ int GModel::readUNV(const std::string &name) ...@@ -1186,24 +1186,77 @@ int GModel::readUNV(const std::string &name)
if(!fscanf(fp, "%d", &n[i])) return 0; if(!fscanf(fp, "%d", &n[i])) return 0;
if(!checkVertexIndex(n[i], vertices)) return 0; if(!checkVertexIndex(n[i], vertices)) return 0;
} }
int type_msh = 0;
int dim = -1;
switch(type){ switch(type){
case 11: case 21: case 22: case 31: type_msh = LGN1; break; case 11: case 21: case 22: case 31:
case 23: case 24: case 32: type_msh = LGN2; break; elements[0][elementary].push_back
case 41: case 51: case 61: case 74: case 81: case 91: type_msh = TRI1; break; (new MLine(vertices[n[0]], vertices[n[1]], num));
case 42: case 52: case 62: case 72: case 82: case 92: type_msh = TRI2; break; dim = 1;
case 44: case 54: case 64: case 71: case 84: case 94: type_msh = QUA1; break; break;
case 45: case 55: case 65: case 75: case 85: case 95: type_msh = QUA2; break; case 23: case 24: case 32:
case 111: type_msh = TET1; break; elements[0][elementary].push_back
case 118: type_msh = TET2; break; (new MLine3(vertices[n[0]], vertices[n[2]], vertices[n[1]], num));
case 104: case 115: type_msh = HEX1; break; dim = 1;
case 105: case 116: type_msh = HEX2; break; break;
case 101: case 112: type_msh = PRI1; break; case 41: case 51: case 61: case 74: case 81: case 91:
case 102: case 113: type_msh = PRI2; break; elements[1][elementary].push_back
(new MTriangle(vertices[n[0]], vertices[n[1]], vertices[n[2]], num));
dim = 2;
break;
case 42: case 52: case 62: case 72: case 82: case 92:
elements[1][elementary].push_back
(new MTriangle6(vertices[n[0]], vertices[n[2]], vertices[n[4]],
vertices[n[1]], vertices[n[3]], vertices[n[5]], num));
dim = 2;
break;
case 44: case 54: case 64: case 71: case 84: case 94:
elements[2][elementary].push_back
(new MQuadrangle(vertices[n[0]], vertices[n[1]], vertices[n[2]],
vertices[n[3]], num));
dim = 2;
break;
case 45: case 55: case 65: case 75: case 85: case 95:
Msg(WARNING, "Quad8 is not implemented yet");
break;
case 111:
elements[3][elementary].push_back
(new MTetrahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]],
vertices[n[3]], num));
dim = 3;
break;
case 118:
elements[3][elementary].push_back
(new MTetrahedron10(vertices[n[0]], vertices[n[2]], vertices[n[4]],
vertices[n[9]], vertices[n[1]], vertices[n[3]],
vertices[n[5]], vertices[n[8]], vertices[n[6]],
vertices[n[7]], num));
dim = 3;
break;
case 104: case 115:
elements[4][elementary].push_back
(new MHexahedron(vertices[n[0]], vertices[n[1]], vertices[n[2]],
vertices[n[3]], vertices[n[4]], vertices[n[5]],
vertices[n[6]], vertices[n[7]], num));
dim = 3;
break;
case 105: case 116:
Msg(WARNING, "Hexa20 is not implemented yet");
break;
case 101: case 112:
elements[5][elementary].push_back
(new MPrism(vertices[n[0]], vertices[n[1]], vertices[n[2]],
vertices[n[3]], vertices[n[4]], vertices[n[5]], num));
dim = 3;
break;
case 102: case 113:
Msg(WARNING, "Prism15 is not implemented yet");
break;
} }
if(type_msh && getNumVerticesForElementTypeMSH(type_msh) == numNodes)
createElementMSH(this, num, type_msh, physical, elementary, 0, if(dim >= 0 && physical && (!physicals[dim].count(elementary) ||
n, vertices, points, elements, physicals); !physicals[dim][elementary].count(physical)))
physicals[dim][elementary][physical] = "unnamed";
} }
} }
} }
......
Geo/MElement.cpp
+ 2
− 2
View file @ c3be78b0
// $Id: MElement.cpp,v 1.15 2006-09-03 07:44:10 geuzaine Exp $ // $Id: MElement.cpp,v 1.16 2006-09-07 19:45:15 geuzaine Exp $
// //
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
// //
...@@ -264,7 +264,7 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical) ...@@ -264,7 +264,7 @@ void MElement::writeUNV(FILE *fp, int num, int elementary, int physical)
if(type == 21 || type == 24) // linear beam or parabolic beam if(type == 21 || type == 24) // linear beam or parabolic beam
fprintf(fp, "%10d%10d%10d\n", 0, 0, 0); fprintf(fp, "%10d%10d%10d\n", 0, 0, 0);
for(int k = 0; k < n; k++) { for(int k = 0; k < n; k++) {
fprintf(fp, "%10d", getVertex(k)->getNum()); fprintf(fp, "%10d", getVertexUNV(k)->getNum());
if(k % 8 == 7) if(k % 8 == 7)
fprintf(fp, "\n"); fprintf(fp, "\n");
} }
......
Geo/MElement.h
+ 41
− 11
View file @ c3be78b0
...@@ -89,6 +89,9 @@ class MElement ...@@ -89,6 +89,9 @@ class MElement
virtual int getNumVertices() = 0; virtual int getNumVertices() = 0;
virtual MVertex *getVertex(int num) = 0; virtual MVertex *getVertex(int num) = 0;
// get the vertex using the I-deas UNV ordering
virtual MVertex *getVertexUNV(int num){ return getVertex(num); }
// get the number of vertices associated with edges, faces and // get the number of vertices associated with edges, faces and
// volumes (nonzero only for higher order elements) // volumes (nonzero only for higher order elements)
virtual int getNumEdgeVertices(){ return 0; } virtual int getNumEdgeVertices(){ return 0; }
...@@ -172,29 +175,34 @@ class MLine : public MElement { ...@@ -172,29 +175,34 @@ class MLine : public MElement {
virtual char *getStringForBDF(){ return "CBAR"; } virtual char *getStringForBDF(){ return "CBAR"; }
}; };
class MLine2 : public MLine { class MLine3 : public MLine {
protected: protected:
MVertex *_vs[1]; MVertex *_vs[1];
public : public :
MLine2(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0) MLine3(MVertex *v0, MVertex *v1, MVertex *v2, int num=0, int part=0)
: MLine(v0, v1, num, part) : MLine(v0, v1, num, part)
{ {
_vs[0] = v2; _vs[0] = v2;
} }
MLine2(std::vector<MVertex*> &v, int num=0, int part=0) MLine3(std::vector<MVertex*> &v, int num=0, int part=0)
: MLine(v, num, part) : MLine(v, num, part)
{ {
_vs[0] = v[2]; _vs[0] = v[2];
} }
~MLine2(){} ~MLine3(){}
virtual int getPolynomialOrder(){ return 2; } virtual int getPolynomialOrder(){ return 2; }
virtual int getNumVertices(){ return 3; } virtual int getNumVertices(){ return 3; }
virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; } virtual MVertex *getVertex(int num){ return num < 2 ? _v[num] : _vs[num - 2]; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[3] = {0, 2, 1};
return getVertex(map[num]);
}
virtual int getNumEdgeVertices(){ return 1; } virtual int getNumEdgeVertices(){ return 1; }
virtual int getNumEdgesRep(){ return 2; } virtual int getNumEdgesRep(){ return 2; }
virtual MEdge getEdgeRep(int num) virtual MEdge getEdgeRep(int num)
{ {
static int edges_lin2[2][2] = { static const int edges_lin2[2][2] = {
{0, 2}, {2, 1}, {0, 2}, {2, 1},
}; };
int i0 = edges_lin2[num][0]; int i0 = edges_lin2[num][0];
...@@ -266,11 +274,16 @@ class MTriangle6 : public MTriangle { ...@@ -266,11 +274,16 @@ class MTriangle6 : public MTriangle {
virtual int getPolynomialOrder(){ return 2; } virtual int getPolynomialOrder(){ return 2; }
virtual int getNumVertices(){ return 6; } virtual int getNumVertices(){ return 6; }
virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; } virtual MVertex *getVertex(int num){ return num < 3 ? _v[num] : _vs[num - 3]; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[6] = {0, 3, 1, 4, 2, 5};
return getVertex(map[num]);
}
virtual int getNumEdgeVertices(){ return 3; } virtual int getNumEdgeVertices(){ return 3; }
virtual int getNumEdgesRep(){ return 6; } virtual int getNumEdgesRep(){ return 6; }
virtual MEdge getEdgeRep(int num) virtual MEdge getEdgeRep(int num)
{ {
static int edges_tri2[6][2] = { static const int edges_tri2[6][2] = {
{0, 3}, {3, 1}, {0, 3}, {3, 1},
{1, 4}, {4, 2}, {1, 4}, {4, 2},
{2, 5}, {5, 0}, {2, 5}, {5, 0},
...@@ -282,7 +295,7 @@ class MTriangle6 : public MTriangle { ...@@ -282,7 +295,7 @@ class MTriangle6 : public MTriangle {
virtual int getNumFacesRep(){ return 4; } virtual int getNumFacesRep(){ return 4; }
virtual MFace getFaceRep(int num) virtual MFace getFaceRep(int num)
{ {
static int trifaces_tri2[4][3] = { static const int trifaces_tri2[4][3] = {
{0, 3, 5}, {0, 3, 5},
{1, 4, 3}, {1, 4, 3},
{2, 5, 4}, {2, 5, 4},
...@@ -355,11 +368,15 @@ class MQuadrangle9 : public MQuadrangle { ...@@ -355,11 +368,15 @@ class MQuadrangle9 : public MQuadrangle {
virtual int getNumFaceVertices(){ return 1; } virtual int getNumFaceVertices(){ return 1; }
// TODO: edgeRep, faceRep // TODO: edgeRep, faceRep
virtual int getTypeForMSH(){ return QUA2; } virtual int getTypeForMSH(){ return QUA2; }
virtual int getTypeForUNV(){ return 95; } // thin shell parabolic quadrilateral virtual int getTypeForUNV(){ return 0; } // not available
virtual char *getStringForPOS(){ return "SQ2"; } virtual char *getStringForPOS(){ return "SQ2"; }
virtual char *getStringForBDF(){ return 0; } // not available virtual char *getStringForBDF(){ return 0; } // not available
}; };
// TODO: MQuadrangle8
// virtual int getTypeForUNV(){ return 95; } // shell parabolic quadrilateral
// virtual char *getStringForBDF(){ return "CQUAD8"; }
class MTetrahedron : public MElement { class MTetrahedron : public MElement {
protected: protected:
MVertex *_v[4]; MVertex *_v[4];
...@@ -440,12 +457,17 @@ class MTetrahedron10 : public MTetrahedron { ...@@ -440,12 +457,17 @@ class MTetrahedron10 : public MTetrahedron {
virtual int getPolynomialOrder(){ return 2; } virtual int getPolynomialOrder(){ return 2; }
virtual int getNumVertices(){ return 10; } virtual int getNumVertices(){ return 10; }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; } virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
virtual MVertex *getVertexUNV(int num)
{
static const int map[10] = {0, 4, 1, 5, 2, 6, 8, 9, 7, 3};
return getVertex(map[num]);
}
virtual int getNumEdgeVertices(){ return 6; } virtual int getNumEdgeVertices(){ return 6; }
// TODO: edgeRep, faceRep // TODO: edgeRep, faceRep
virtual int getTypeForMSH(){ return TET2; } virtual int getTypeForMSH(){ return TET2; }
virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
virtual char *getStringForPOS(){ return "SS2"; } virtual char *getStringForPOS(){ return "SS2"; }
virtual char *getStringForBDF(){ return 0; } // not available virtual char *getStringForBDF(){ return 0; } // TODO: "CTETRA"
virtual void setVolumePositive() virtual void setVolumePositive()
{ {
if(getVolumeSign() < 0){ if(getVolumeSign() < 0){
...@@ -549,7 +571,7 @@ class MHexahedron27 : public MHexahedron { ...@@ -549,7 +571,7 @@ class MHexahedron27 : public MHexahedron {
virtual int getNumVolumeVertices(){ return 1; } virtual int getNumVolumeVertices(){ return 1; }
// TODO: edgeRep, faceRep // TODO: edgeRep, faceRep
virtual int getTypeForMSH(){ return HEX2; } virtual int getTypeForMSH(){ return HEX2; }
virtual int getTypeForUNV(){ return 116; } // solid parabolic brick virtual int getTypeForUNV(){ return 0; } // not available
virtual char *getStringForPOS(){ return "SH2"; } virtual char *getStringForPOS(){ return "SH2"; }
virtual char *getStringForBDF(){ return 0; } // not available virtual char *getStringForBDF(){ return 0; } // not available
virtual void setVolumePositive() virtual void setVolumePositive()
...@@ -568,6 +590,10 @@ class MHexahedron27 : public MHexahedron { ...@@ -568,6 +590,10 @@ class MHexahedron27 : public MHexahedron {
} }
}; };
// TODO: MHexahedron20
// virtual int getTypeForUNV(){ return 116; } // solid parabolic brick
// virtual char *getStringForBDF(){ return "CHEXA"; }
class MPrism : public MElement { class MPrism : public MElement {
protected: protected:
MVertex *_v[6]; MVertex *_v[6];
...@@ -661,7 +687,7 @@ class MPrism18 : public MPrism { ...@@ -661,7 +687,7 @@ class MPrism18 : public MPrism {
virtual int getNumFaceVertices(){ return 3; } virtual int getNumFaceVertices(){ return 3; }
// TODO: edgeRep, faceRep // TODO: edgeRep, faceRep
virtual int getTypeForMSH(){ return PRI2; } virtual int getTypeForMSH(){ return PRI2; }
virtual int getTypeForUNV(){ return 113; } // solid parabolic wedge virtual int getTypeForUNV(){ return 0; } // not available
virtual char *getStringForPOS(){ return "SI2"; } virtual char *getStringForPOS(){ return "SI2"; }
virtual char *getStringForBDF(){ return 0; } // not available virtual char *getStringForBDF(){ return 0; } // not available
virtual void setVolumePositive() virtual void setVolumePositive()
...@@ -677,6 +703,10 @@ class MPrism18 : public MPrism { ...@@ -677,6 +703,10 @@ class MPrism18 : public MPrism {
} }
}; };
// TODO: MPrism15
// virtual int getTypeForUNV(){ return 113; } // solid parabolic wedge
// virtual char *getStringForBDF(){ return "CPENTA"; }
class MPyramid : public MElement { class MPyramid : public MElement {
protected: protected:
MVertex *_v[5]; MVertex *_v[5];
......
Geo/fourierModel.cpp
+ 40
− 24
View file @ c3be78b0
...@@ -38,36 +38,29 @@ public: ...@@ -38,36 +38,29 @@ public:
} }
}; };
static void getParamForPointInOverlaps(GModel *m, std::vector<int> &overlaps,
double x, double y, double z,
std::vector<std::pair<GFace*, SPoint2> > &param)
{
// we only loop on patches that are known a priori to overlap with
// the current patch (this way we don't normalize the POU across
// hard edges)
const double tol = 1.e-2; // need to adapt this w.r.t size of model
for(unsigned int i = 0; i < overlaps.size(); i++){
GFace *f = m->faceByTag(overlaps[i]);
SPoint2 uv = f->parFromPoint(SPoint3(x, y, z));
if(f->containsParam(uv)){
GPoint xyz = f->point(uv);
if(fabs(xyz.x() - x) < tol && fabs(xyz.y() - y) < tol && fabs(xyz.z() - z) < tol)
param.push_back(std::pair<GFace*, SPoint2>(f, uv));
}
}
}
class computePartitionOfUnity{ class computePartitionOfUnity{
public: public:
void operator() (GFace *gf) void operator() (GFace *gf)
{ {
// we only normalize the partition of unity amongst patches that
// overlap; we don't normalize across hard edges
std::vector<int> overlaps; std::vector<int> overlaps;
FM->GetNeighbor(gf->tag(), overlaps); FM->GetNeighbor(gf->tag(), overlaps);
for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++){ for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
MVertex *v = gf->mesh_vertices[i]; MVertex *v = gf->mesh_vertices[i];
std::vector<std::pair<GFace*, SPoint2> > param; std::vector<std::pair<GFace*, SPoint2> > param;
getParamForPointInOverlaps(gf->model(), overlaps, v->x(), v->y(), v->z(), param); for(unsigned int j = 0; j < overlaps.size(); j++){
GFace *gf2 = gf->model()->faceByTag(overlaps[j]);
SPoint2 uv2 = gf2->parFromPoint(SPoint3(v->x(), v->y(), v->z()));
if(gf2->containsParam(uv2)){
GPoint xyz2 = gf2->point(uv2);
const double tol = 1.e-2; // need to adapt this w.r.t size of model
if(fabs(xyz2.x() - v->x()) < tol &&
fabs(xyz2.y() - v->y()) < tol &&
fabs(xyz2.z() - v->z()) < tol)
param.push_back(std::pair<GFace*, SPoint2>(gf2, uv2));
}
}
if(param.empty()){ if(param.empty()){
Msg(GERROR, "Vertex %d does not belong to any patch", v->getNum()); Msg(GERROR, "Vertex %d does not belong to any patch", v->getNum());
} }
...@@ -91,7 +84,7 @@ public: ...@@ -91,7 +84,7 @@ public:
} }
}; };
class removeGrout{ class createGrooves{
public: public:
void operator() (GFace *gf) void operator() (GFace *gf)
{ {
...@@ -108,6 +101,27 @@ public: ...@@ -108,6 +101,27 @@ public:
} }
}; };
class createGrout{
public:
void operator() (GFace *gf)
{
printf("processing grout for face %d\n", gf->tag());
std::vector<int> overlaps;
FM->GetNeighbor(gf->tag(), overlaps);
for(unsigned int j = 0; j < overlaps.size(); j++){
GFace *gf2 = gf->model()->faceByTag(overlaps[j]);
if(gf != gf2){
}
}
}
};
class exportFourierFace{ class exportFourierFace{
public: public:
void operator() (GFace *gf) void operator() (GFace *gf)
...@@ -155,9 +169,11 @@ fourierModel::fourierModel(const std::string &name) ...@@ -155,9 +169,11 @@ fourierModel::fourierModel(const std::string &name)
// compute partition of unity // compute partition of unity
std::for_each(firstFace(), lastFace(), computePartitionOfUnity()); std::for_each(firstFace(), lastFace(), computePartitionOfUnity());
// remove the grout // create grooves
std::for_each(firstFace(), lastFace(), removeGrout()); std::for_each(firstFace(), lastFace(), createGrooves());
// create grout
std::for_each(firstFace(), lastFace(), createGrout());
// visualize as a post-pro view // visualize as a post-pro view
std::for_each(firstFace(), lastFace(), exportFourierFace()); std::for_each(firstFace(), lastFace(), exportFourierFace());
......
Mesh/SecondOrder.cpp
+ 2
− 2
View file @ c3be78b0
// $Id: SecondOrder.cpp,v 1.44 2006-09-04 15:58:22 geuzaine Exp $ // $Id: SecondOrder.cpp,v 1.45 2006-09-07 19:45:15 geuzaine Exp $
// //
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle // Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
// //
...@@ -195,7 +195,7 @@ void setSecondOrder(GEdge *ge, ...@@ -195,7 +195,7 @@ void setSecondOrder(GEdge *ge,
MLine *l = ge->lines[i]; MLine *l = ge->lines[i];
std::vector<MVertex*> ve; std::vector<MVertex*> ve;
getEdgeVertices(ge, l, ve, edgeVertices, linear); getEdgeVertices(ge, l, ve, edgeVertices, linear);
lines2.push_back(new MLine2(l->getVertex(0), l->getVertex(1), ve[0])); lines2.push_back(new MLine3(l->getVertex(0), l->getVertex(1), ve[0]));
delete l; delete l;
} }
ge->lines = lines2; ge->lines = lines2;
......
doc/CREDITS
+ 2
− 2
View file @ c3be78b0
$Id: CREDITS,v 1.36 2006-08-10 15:29:26 geuzaine Exp $ $Id: CREDITS,v 1.37 2006-09-07 19:45:15 geuzaine Exp $
Gmsh is copyright (C) 1997-2006 Gmsh is copyright (C) 1997-2006
...@@ -117,4 +117,4 @@ Krishna Mohan Gundu <gkmohan at gmail.com>, Christopher Stott <C.Stott ...@@ -117,4 +117,4 @@ Krishna Mohan Gundu <gkmohan at gmail.com>, Christopher Stott <C.Stott
at surrey.ac.uk>, Timmy Schumacher <Tim.Schumacher at colorado.edu>, at surrey.ac.uk>, Timmy Schumacher <Tim.Schumacher at colorado.edu>,
Carl Osterwisch <osterwischc at asme.org>, Bruno Frackowiak Carl Osterwisch <osterwischc at asme.org>, Bruno Frackowiak
<bruno.frackowiak at onecert.fr>, Philip Kelleners <P.H.Kelleners at <bruno.frackowiak at onecert.fr>, Philip Kelleners <P.H.Kelleners at
ctw.utwente.nl>. ctw.utwente.nl>, Romuald Conty <conty@phimeca.com>.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment