diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index d78f5aa7f2873affa7882bc5b0bf2088a4debe03..8afa3f41bb895b14d4d5172c482ee15084b78ef6 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -82,6 +82,8 @@
#define MSH_TET_20 29
#define MSH_TET_35 30
#define MSH_TET_56 31
+#define MSH_TET_34 32
+#define MSH_TET_52 33
// Geometric entities
#define ENT_NONE 0
diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 6cf735bea49dcbc8790d1afb0c6f79dd8436037e..6f79f03ff0b5e437a33e5fe7f7037999791ad732 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.69 2008-07-10 13:29:24 geuzaine Exp $
+// $Id: GFace.cpp,v 1.70 2008-07-11 10:54:24 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -732,7 +732,7 @@ bool GFace::buildSTLTriangulation()
// by default we assume that straight lines are geodesics
SPoint2 GFace::geodesic(const SPoint2 &pt1 , const SPoint2 &pt2 , double t)
{
- if(CTX.mesh.second_order_experimental){
+ if(CTX.mesh.second_order_experimental && geomType() != GEntity::Plane){
// FIXME: this is buggy -- remove the CTX option once we do it in
// a robust manner
GPoint gp1 = point(pt1.x(), pt1.y());
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index e13d701e10514933d18d6de7c906f4b1ec7c1b56..da743dac7a2ffced3acfb0c4f34a0aa20d86c24b 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.58 2008-07-08 12:44:33 remacle Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.59 2008-07-11 10:54:24 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -140,8 +140,9 @@ static int getNumVerticesForElementTypeMSH(int type)
case MSH_TET_4 : return 4;
case MSH_TET_10 : return 4 + 6;
case MSH_TET_20 : return 4 + 12 + 4;
- case MSH_TET_35 : return 4 + 16 + 12 + 3;
- case MSH_TET_56 : return 4 + 20 + 24 + 8;
+ case MSH_TET_35 : return 4 + 18 + 12 + 1;
+ case MSH_TET_52 : return 4 + 24 + 24 + 0;
+ case MSH_TET_56 : return 4 + 24 + 24 + 4;
case MSH_HEX_8 : return 8;
case MSH_HEX_20 : return 8 + 12;
case MSH_HEX_27 : return 8 + 12 + 6 + 1;
@@ -171,6 +172,7 @@ static void createElementMSH(GModel *m, int num, int type, int physical,
else{
MElementFactory factory;
MElement *e = factory.create(type, v, num, part);
+
if(!e){
Msg::Error("Unknown type of element %d", type);
return;
@@ -333,7 +335,7 @@ int GModel::readMSH(const std::string &name)
}
if(!(numVertices = getNumVerticesForElementTypeMSH(type))) return 0;
}
- int indices[30];
+ int indices[60];
for(int j = 0; j < numVertices; j++) fscanf(fp, "%d", &indices[j]);
std::vector<MVertex*> vertices;
if(vertexVector.size()){
@@ -447,7 +449,7 @@ int GModel::readMSH(const std::string &name)
storeVerticesInEntities(vertexMap);
// store the physical tags
- for(int i = 0; i < 4; i++)
+ for(int i = 0; i < 4; i++)
storePhysicalTagsInEntities(this, i, physicals[i]);
fclose(fp);
@@ -634,7 +636,7 @@ int GModel::writeMSH(const std::string &name, double version, bool binary,
for(fiter it = firstFace(); it != lastFace(); ++it)
writeElementsMSH(fp, (*it)->quadrangles, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
- writeElementHeaderMSH(binary, fp, elements, MSH_TET_4, MSH_TET_10);
+ writeElementHeaderMSH(binary, fp, elements, MSH_TET_4, MSH_TET_10, MSH_TET_20, MSH_TET_35, MSH_TET_56, MSH_TET_52);
for(riter it = firstRegion(); it != lastRegion(); ++it)
writeElementsMSH(fp, (*it)->tetrahedra, saveAll, version, binary, num,
(*it)->tag(), (*it)->physicals);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6a373763a5986a3378ce0d65b47b3146d7afdf23..a8168ffe6ae2be730bde8826281d849c518606fb 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.79 2008-07-10 13:29:24 geuzaine Exp $
+// $Id: MElement.cpp,v 1.80 2008-07-11 10:54:24 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -418,7 +418,7 @@ void MElement::writeMSH(FILE *fp, double version, bool binary, int num,
if(physical < 0) revert();
- int verts[30];
+ int verts[60];
for(int i = 0; i < n; i++)
verts[i] = getVertex(i)->getIndex();
@@ -553,7 +553,7 @@ void MElement::writeVTK(FILE *fp, bool binary)
int n = getNumVertices();
if(binary){
- int verts[30];
+ int verts[60];
verts[0] = n;
for(int i = 0; i < n; i++)
verts[i + 1] = getVertexVTK(i)->getIndex() - 1;
@@ -731,13 +731,25 @@ void MTetrahedron::pnt(int ord, MVertex *vs[], double uu, double vv, double ww,S
{
#if !defined(HAVE_GMSH_EMBEDDED)
double sf[256];
- switch(ord){
- case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, ww,sf); break;
- case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww,sf); break;
- case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww,sf); break;
- case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww,sf); break;
- case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww,sf); break;
- default: Msg::Error("Order %d tetrahedron pnt not implemented", ord); break;
+
+ int nv = getNumVolumeVertices();
+
+ if (!nv){
+ switch(ord){
+ case 1: gmshFunctionSpaces::find(MSH_TET_4).f(uu, vv, ww,sf); break;
+ case 2: gmshFunctionSpaces::find(MSH_TET_10).f(uu, vv, ww,sf); break;
+ case 3: gmshFunctionSpaces::find(MSH_TET_20).f(uu, vv, ww,sf); break;
+ case 4: gmshFunctionSpaces::find(MSH_TET_34).f(uu, vv, ww,sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_52).f(uu, vv, ww,sf); break;
+ default: Msg::Error("Order %d tetrahedron pnt not implemented", ord); break;
+ }
+ }
+ else{
+ switch(ord){
+ case 4: gmshFunctionSpaces::find(MSH_TET_35).f(uu, vv, ww,sf); break;
+ case 5: gmshFunctionSpaces::find(MSH_TET_56).f(uu, vv, ww,sf); break;
+ default: Msg::Error("Order %d tetrahedron pnt not implemented", ord); break;
+ }
}
double x = 0 ; for(int i = 0; i < 4; i++) x += sf[i] * _v[i]->x();
@@ -884,7 +896,7 @@ void MTetrahedron::jac( double uu, double vv, double ww, double j[3][3]) {
return jac(1,0,uu,vv,ww,j);
}
-const int numSubEdges = 12;
+const int numSubEdges = 6;
int MTriangleN::getNumFacesRep(){ return numSubEdges * numSubEdges; }
@@ -1020,7 +1032,9 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_PYR_13: return new MPyramid13(v, num, part);
case MSH_PYR_14: return new MPyramid14(v, num, part);
case MSH_TET_20: return new MTetrahedronN(v, 3, num, part);
+ case MSH_TET_34: return new MTetrahedronN(v, 3, num, part);
case MSH_TET_35: return new MTetrahedronN(v, 4, num, part);
+ case MSH_TET_52: return new MTetrahedronN(v, 5, num, part);
case MSH_TET_56: return new MTetrahedronN(v, 5, num, part);
default: return 0;
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index cf4e67e4724aefc841e4e409ce5c82b18f2e9228..ea226c1b6cb41fb7aef293f5293c4f3583a20e82 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -1479,12 +1479,24 @@ class MTetrahedronN : public MTetrahedron {
const int ie = (num+1)*getNumFaceVertices();
for(int i = num*getNumFaceVertices(); i != ie; ++i) v[j++] = _vs[i];
}
+ virtual int getNumVolumeVertices()
+ {
+ switch(getTypeForMSH()){
+ case MSH_TET_35 : return 1;
+ case MSH_TET_56 : return 4;
+ default : return 0;
+ }
+ }
+ virtual int getNumEdgeVertices(){ return _order - 1; }
+
virtual int getTypeForMSH()
{
// (p+1)*(p+2)*(p+3)/6
if(_order == 3 && _vs.size() + 4 == 20) return MSH_TET_20;
+ if(_order == 4 && _vs.size() + 4 == 34) return MSH_TET_34;
if(_order == 4 && _vs.size() + 4 == 35) return MSH_TET_35;
if(_order == 5 && _vs.size() + 4 == 56) return MSH_TET_56;
+ if(_order == 5 && _vs.size() + 4 == 52) return MSH_TET_52;
return 0;
}
virtual void revert()
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 5d30f79b006cc2db7b2a467821faa62d62b1c515..c71d68d7a7c4a33999b0e43f91260407f178fc76 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1,4 +1,4 @@
-// $Id: HighOrder.cpp,v 1.34 2008-07-10 13:29:24 geuzaine Exp $
+// $Id: HighOrder.cpp,v 1.35 2008-07-11 10:54:24 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -650,6 +650,51 @@ void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
}
}
+
+void getRegionVertices(GRegion *gr,
+ MElement *incomplete,
+ MElement *ele,
+ std::vector<MVertex*> &vr,
+ bool linear, int nPts = 1)
+{
+
+ Double_Matrix points;
+ int start = 0;
+
+ switch (nPts){
+ case 3 :
+ points = gmshFunctionSpaces::find(MSH_TET_35).points;
+ start = 34;
+ break;
+ case 4 :
+ points = gmshFunctionSpaces::find(MSH_TET_56).points;
+ start = 52;
+ break;
+ default :
+ return;
+ break;
+ }
+
+ for(int k = start ; k < points.size1() ; k++){
+ MVertex *v;
+ const double t1 = points(k, 0);
+ const double t2 = points(k, 1);
+ const double t3 = points(k, 2);
+ double X(0),Y(0),Z(0);
+ for (int j=0; j<incomplete->getNumVertices(); j++){
+ double sf ; incomplete->getShapeFunction(j,t1,t2,t3,sf);
+ MVertex *vt = incomplete->getVertex(j);
+ X += sf * vt->x();
+ Y += sf * vt->y();
+ Z += sf * vt->z();
+ }
+ v = new MVertex(X,Y,Z, gr);
+ gr->mesh_vertices.push_back(v);
+ vr.push_back(v);
+ }
+}
+
+
void setHighOrder(GEdge *ge, edgeContainer &edgeVertices, bool linear,
int nbPts = 1)
{
@@ -736,7 +781,7 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
std::vector<MTetrahedron*> tetrahedra2;
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
MTetrahedron *t = gr->tetrahedra[i];
- std::vector<MVertex*> ve,vf;
+ std::vector<MVertex*> ve,vf,vr;
getEdgeVertices(gr, t, ve, edgeVertices, linear, nPts);
if (nPts == 1){
tetrahedra2.push_back
@@ -745,7 +790,11 @@ void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
}
else{
getFaceVertices(gr, t, vf, faceVertices, linear, nPts);
- ve.insert(ve.end(), vf.begin(), vf.end());
+ ve.insert(ve.end(), vf.begin(), vf.end());
+ MTetrahedronN incpl (t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
+ ve, nPts + 1);
+ getRegionVertices(gr, &incpl, t, vr, linear, nPts);
+ ve.insert(ve.end(), vr.begin(), vr.end());
tetrahedra2.push_back
(new MTetrahedronN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
ve, nPts + 1));
@@ -1432,7 +1481,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete)
int nPts = order - 1;
- Msg::StatusBar(1, true, "Meshing second order...");
+ Msg::StatusBar(1, true, "Generating High Order Nodes (q = %d) ...",order);
double t1 = Cpu();
// first, make sure to remove any existsing second order vertices/elements
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index bebaa356a197e654f3067912875f663d8182f0b8..71646f670420c662a1632b0908715088270b6b7e 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -1,4 +1,4 @@
-// $Id: FunctionSpace.cpp,v 1.8 2008-07-03 17:06:04 geuzaine Exp $
+// $Id: FunctionSpace.cpp,v 1.9 2008-07-11 10:54:24 remacle Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -613,6 +613,14 @@ const gmshFunctionSpace &gmshFunctionSpaces::find(int tag)
F.monomials = generatePascalTetrahedron(4);
F.points = gmshGeneratePointsTetrahedron(4, false);
break;
+ case MSH_TET_34 :
+ F.monomials = generatePascalSerendipityTetrahedron(4);
+ F.points = gmshGeneratePointsTetrahedron(4, true);
+ break;
+ case MSH_TET_52 :
+ F.monomials = generatePascalSerendipityTetrahedron(5);
+ F.points = gmshGeneratePointsTetrahedron(5, true);
+ break;
case MSH_TET_56 :
F.monomials = generatePascalTetrahedron(5);
F.points = gmshGeneratePointsTetrahedron(5, false);