From c3be78b09edd907165e3df4a89abc2842fc419ad Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Thu, 7 Sep 2006 19:45:15 +0000
Subject: [PATCH] - Use correct vertex ordering for second order simplices in
UNV format (thanks to Romuald Conty <conty@phimeca.com>)
- MLine2 -> MLine3
---
Geo/GModelIO.cpp | 89 +++++++++++++++++++++++++++++++++++---------
Geo/MElement.cpp | 4 +-
Geo/MElement.h | 52 ++++++++++++++++++++------
Geo/fourierModel.cpp | 64 +++++++++++++++++++------------
Mesh/SecondOrder.cpp | 4 +-
doc/CREDITS | 4 +-
6 files changed, 158 insertions(+), 59 deletions(-)
diff --git a/Geo/GModelIO.cpp b/Geo/GModelIO.cpp
index bb1d95e189..a49769ed0d 100644
--- a/Geo/GModelIO.cpp
+++ b/Geo/GModelIO.cpp
@@ -1,4 +1,4 @@
-// $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
//
@@ -236,7 +236,7 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
break;
case LGN2:
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;
break;
case TRI1:
@@ -1186,24 +1186,77 @@ int GModel::readUNV(const std::string &name)
if(!fscanf(fp, "%d", &n[i])) return 0;
if(!checkVertexIndex(n[i], vertices)) return 0;
}
- int type_msh = 0;
+
+ int dim = -1;
switch(type){
- case 11: case 21: case 22: case 31: type_msh = LGN1; break;
- case 23: case 24: case 32: type_msh = LGN2; break;
- case 41: case 51: case 61: case 74: case 81: case 91: type_msh = TRI1; break;
- case 42: case 52: case 62: case 72: case 82: case 92: type_msh = TRI2; break;
- case 44: case 54: case 64: case 71: case 84: case 94: type_msh = QUA1; break;
- case 45: case 55: case 65: case 75: case 85: case 95: type_msh = QUA2; break;
- case 111: type_msh = TET1; break;
- case 118: type_msh = TET2; break;
- case 104: case 115: type_msh = HEX1; break;
- case 105: case 116: type_msh = HEX2; break;
- case 101: case 112: type_msh = PRI1; break;
- case 102: case 113: type_msh = PRI2; break;
+ case 11: case 21: case 22: case 31:
+ elements[0][elementary].push_back
+ (new MLine(vertices[n[0]], vertices[n[1]], num));
+ dim = 1;
+ break;
+ case 23: case 24: case 32:
+ elements[0][elementary].push_back
+ (new MLine3(vertices[n[0]], vertices[n[2]], vertices[n[1]], num));
+ dim = 1;
+ break;
+ case 41: case 51: case 61: case 74: case 81: case 91:
+ 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,
- n, vertices, points, elements, physicals);
+
+ if(dim >= 0 && physical && (!physicals[dim].count(elementary) ||
+ !physicals[dim][elementary].count(physical)))
+ physicals[dim][elementary][physical] = "unnamed";
}
}
}
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index b143f604be..08ce1e8bff 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $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
//
@@ -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
fprintf(fp, "%10d%10d%10d\n", 0, 0, 0);
for(int k = 0; k < n; k++) {
- fprintf(fp, "%10d", getVertex(k)->getNum());
+ fprintf(fp, "%10d", getVertexUNV(k)->getNum());
if(k % 8 == 7)
fprintf(fp, "\n");
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index daba42afc1..e24ef8e4a8 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -89,6 +89,9 @@ class MElement
virtual int getNumVertices() = 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
// volumes (nonzero only for higher order elements)
virtual int getNumEdgeVertices(){ return 0; }
@@ -172,29 +175,34 @@ class MLine : public MElement {
virtual char *getStringForBDF(){ return "CBAR"; }
};
-class MLine2 : public MLine {
+class MLine3 : public MLine {
protected:
MVertex *_vs[1];
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)
{
_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)
{
_vs[0] = v[2];
}
- ~MLine2(){}
+ ~MLine3(){}
virtual int getPolynomialOrder(){ return 2; }
virtual int getNumVertices(){ return 3; }
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 getNumEdgesRep(){ return 2; }
virtual MEdge getEdgeRep(int num)
{
- static int edges_lin2[2][2] = {
+ static const int edges_lin2[2][2] = {
{0, 2}, {2, 1},
};
int i0 = edges_lin2[num][0];
@@ -266,11 +274,16 @@ class MTriangle6 : public MTriangle {
virtual int getPolynomialOrder(){ return 2; }
virtual int getNumVertices(){ return 6; }
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 getNumEdgesRep(){ return 6; }
virtual MEdge getEdgeRep(int num)
{
- static int edges_tri2[6][2] = {
+ static const int edges_tri2[6][2] = {
{0, 3}, {3, 1},
{1, 4}, {4, 2},
{2, 5}, {5, 0},
@@ -282,7 +295,7 @@ class MTriangle6 : public MTriangle {
virtual int getNumFacesRep(){ return 4; }
virtual MFace getFaceRep(int num)
{
- static int trifaces_tri2[4][3] = {
+ static const int trifaces_tri2[4][3] = {
{0, 3, 5},
{1, 4, 3},
{2, 5, 4},
@@ -355,11 +368,15 @@ class MQuadrangle9 : public MQuadrangle {
virtual int getNumFaceVertices(){ return 1; }
// TODO: edgeRep, faceRep
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 *getStringForBDF(){ return 0; } // not available
};
+// TODO: MQuadrangle8
+// virtual int getTypeForUNV(){ return 95; } // shell parabolic quadrilateral
+// virtual char *getStringForBDF(){ return "CQUAD8"; }
+
class MTetrahedron : public MElement {
protected:
MVertex *_v[4];
@@ -440,12 +457,17 @@ class MTetrahedron10 : public MTetrahedron {
virtual int getPolynomialOrder(){ return 2; }
virtual int getNumVertices(){ return 10; }
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; }
// TODO: edgeRep, faceRep
virtual int getTypeForMSH(){ return TET2; }
virtual int getTypeForUNV(){ return 118; } // solid parabolic tetrahedron
virtual char *getStringForPOS(){ return "SS2"; }
- virtual char *getStringForBDF(){ return 0; } // not available
+ virtual char *getStringForBDF(){ return 0; } // TODO: "CTETRA"
virtual void setVolumePositive()
{
if(getVolumeSign() < 0){
@@ -549,7 +571,7 @@ class MHexahedron27 : public MHexahedron {
virtual int getNumVolumeVertices(){ return 1; }
// TODO: edgeRep, faceRep
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 *getStringForBDF(){ return 0; } // not available
virtual void setVolumePositive()
@@ -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 {
protected:
MVertex *_v[6];
@@ -661,7 +687,7 @@ class MPrism18 : public MPrism {
virtual int getNumFaceVertices(){ return 3; }
// TODO: edgeRep, faceRep
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 *getStringForBDF(){ return 0; } // not available
virtual void setVolumePositive()
@@ -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 {
protected:
MVertex *_v[5];
diff --git a/Geo/fourierModel.cpp b/Geo/fourierModel.cpp
index 8d90236ff6..0f8ec65bb0 100644
--- a/Geo/fourierModel.cpp
+++ b/Geo/fourierModel.cpp
@@ -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> > ¶m)
-{
- // 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{
public:
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;
FM->GetNeighbor(gf->tag(), overlaps);
for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
MVertex *v = gf->mesh_vertices[i];
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()){
Msg(GERROR, "Vertex %d does not belong to any patch", v->getNum());
}
@@ -91,7 +84,7 @@ public:
}
};
-class removeGrout{
+class createGrooves{
public:
void operator() (GFace *gf)
{
@@ -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{
public:
void operator() (GFace *gf)
@@ -155,9 +169,11 @@ fourierModel::fourierModel(const std::string &name)
// compute partition of unity
std::for_each(firstFace(), lastFace(), computePartitionOfUnity());
- // remove the grout
- std::for_each(firstFace(), lastFace(), removeGrout());
+ // create grooves
+ std::for_each(firstFace(), lastFace(), createGrooves());
+ // create grout
+ std::for_each(firstFace(), lastFace(), createGrout());
// visualize as a post-pro view
std::for_each(firstFace(), lastFace(), exportFourierFace());
diff --git a/Mesh/SecondOrder.cpp b/Mesh/SecondOrder.cpp
index ff4316d90f..0850575702 100644
--- a/Mesh/SecondOrder.cpp
+++ b/Mesh/SecondOrder.cpp
@@ -1,4 +1,4 @@
-// $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
//
@@ -195,7 +195,7 @@ void setSecondOrder(GEdge *ge,
MLine *l = ge->lines[i];
std::vector<MVertex*> ve;
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;
}
ge->lines = lines2;
diff --git a/doc/CREDITS b/doc/CREDITS
index 74e23949a3..3d1fb040d3 100644
--- a/doc/CREDITS
+++ b/doc/CREDITS
@@ -1,4 +1,4 @@
-$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
@@ -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>,
Carl Osterwisch <osterwischc at asme.org>, Bruno Frackowiak
<bruno.frackowiak at onecert.fr>, Philip Kelleners <P.H.Kelleners at
-ctw.utwente.nl>.
+ctw.utwente.nl>, Romuald Conty <conty@phimeca.com>.
--
GitLab