diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 85613f7800845f81087b2c4e66357622d1925ba1..18b578588860cc0c56ae236fadee15236003cdf1 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -531,8 +531,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
Unew <= umax && Vnew <= vmax &&
Unew >= umin && Vnew >= vmin){
if (onSurface && err2 > 1.e-4 * CTX.lc)
- Msg::Warning("Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g",
- i, j, err, iter, err2);
+ Msg::Warning("Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g in point (%e,%e,%e) on surface %d",
+ i, j, err, iter, err2,X,Y,Z,tag());
return;
}
}
@@ -556,7 +556,7 @@ SPoint2 GFace::parFromPoint(const SPoint3 &p) const
GPoint GFace::closestPoint(const SPoint3 & queryPoint, const double initialGuess[2]) const
{
- Msg::Error("Closet point not implemented for this type of surface");
+ Msg::Error("Closest point not implemented for this type of surface");
return GPoint(0, 0, 0);
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 4ddd476894bcbe0cb4cc907fc474c9d41b28e427..fc96c6362864157706d2a99318498af66165563f 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -13,6 +13,8 @@
#include "MVertex.h"
#include "MEdge.h"
#include "MFace.h"
+#include "Message.h"
+
struct IntPt{
double pt[3];
@@ -1455,6 +1457,67 @@ class MTetrahedron10 : public MTetrahedron {
* `3
*
*/
+
+
+
+/* tet order 3
+
+ * 2
+ * ,/|`\
+ * ,5 | `6 E = order - 1
+ * ,/ 12 `\ C = 4 + 6*E
+ * ,4 | `7 F = ((order - 1)*(order - 2))/2
+ * ,/ | `\ N = total number of vertices
+ * 0-----9-----'.--8-----1
+ * `\. | ,/ Interior vertex numbers
+ * 10. 13 ,14 for edge 0 <= i <= 5: 4+i*E to 3+(i+1)*E
+ * `\. '. 15 for face 0 <= j <= 3: C+j*F to C-1+(j+1)*F
+ * 11\.|/ in volume : C+4*F to N-1
+ * `3
+ *
+ */
+
+
+
+
+static int reverseTet20[20] = {0,2,1,3, // principal vertices
+ 9,8, // E0 switches with E2
+ 7,6, // E1 inverts direction
+ 5,4, // E2 switches with E0
+ 10,11, // E3 pure w edge > remains the same
+ 14,15, // E4 uw edge swithes with v/w edge E5
+ 12,13, // E5 switches with E4
+ 16, // F0 is uv plane, reverts normal
+ 18, // F1 is uw plane, switches with F2
+ 17, // F2 is vw plane, switches with F1
+ 19}; // F3 is uvw plane, reverts normal
+
+static int reverseTet35[35] = {0,2,1,3, // principal vertices
+
+ 12,11,10, // E0 switches with E2
+ 9,8,7, // E1 inverts direction
+ 6,5,4, // E2 switches with E0
+ 13,14,15, // E3 pure w edge > remains the same
+ 19,20,21, // E4 uw edge swithes with v/w edge E5
+ 16,17,18, // E5 switches with E4
+ 22,24,23, // F0 is uv plane, reverts normal
+ 28,30,29, // F1 is uw plane, switches with F2, orientation is different
+ 25,27,26, // F2 is vw plane, switches with F1
+ 31,33,32, // F3 is uvw plane, reverts normal
+ 34}; // central node remains
+
+static int reverseTet34[34] = {0,2,1,3, // principal vertices
+ 12,11,10, // E0 switches with E2
+ 9,8,7, // E1 inverts direction
+ 6,5,4, // E2 switches with E0
+ 13,14,15, // E3 pure w edge > remains the same
+ 19,20,21, // E4 uw edge swithes with v/w edge E5
+ 16,17,18, // E5 switches with E4
+ 22,24,23, // F0 is uv plane, reverts normal
+ 28,29,30, // F1 is uw plane, switches with F2
+ 25,26,27, // F2 is vw plane, switches with F1
+ 31,33,32};// F3 is uvw plane, reverts normal
+
class MTetrahedronN : public MTetrahedron {
protected:
std::vector<MVertex *> _vs;
@@ -1518,12 +1581,38 @@ class MTetrahedronN : public MTetrahedron {
return 0;
}
virtual void revert()
- {
+ {
MVertex *tmp;
- tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
- std::vector<MVertex*> inv;
- inv.insert(inv.begin(), _vs.rbegin(), _vs.rend());
- _vs = inv;
+ tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
+ switch (getTypeForMSH()) {
+ case MSH_TET_20:
+ {
+ std::vector<MVertex*> inv(16);
+ for (int i=0;i<16;i++) inv[i] = _vs[reverseTet20[i+4]-4];
+ _vs = inv;
+ break;
+ }
+ case MSH_TET_35:
+ {
+ std::vector<MVertex*> inv(31);
+ for (int i=0;i<31;i++) inv[i] = _vs[reverseTet35[i+4]-4];
+ _vs = inv;
+ break;
+ }
+ case MSH_TET_34:
+ {
+ std::vector<MVertex*> inv(30);
+ for (int i=0;i<31;i++) inv[i] = _vs[reverseTet34[i+4]-4];
+ _vs = inv;
+ break;
+ }
+ default:
+ {
+ Msg::Error("Reversion of %d order tetrahedron (type %d) not implemented\n",
+ _order,getTypeForMSH());
+ break;
+ }
+ }
}
virtual void jac(double u, double v, double w, double j[3][3])
{
diff --git a/Geo/MFace.cpp b/Geo/MFace.cpp
index ef4520f47ef6cdd1e5ae5496f02886a88a9157b8..3b254cbda380f125c24322ac5451a7eaebcc9aad 100644
--- a/Geo/MFace.cpp
+++ b/Geo/MFace.cpp
@@ -86,3 +86,23 @@ SVector3 MFace::normal() const
return SVector3(n[0], n[1], n[2]);
}
+bool MFace::computeCorrespondence(const MFace& other,int& rotation,bool& swap) const {
+
+ rotation = 0;
+ swap = false;
+
+ if (*this == other) {
+ for (int i=0;i<getNumVertices();i++) {
+ if (_v[0] == other.getVertex(i)) {
+ rotation = i;
+ break;
+ }
+ }
+ if (_v[1] == other.getVertex((rotation+1)%getNumVertices())) swap = false;
+ else swap = true;
+ return true;
+ }
+ return false;
+}
+
+
diff --git a/Geo/MFace.h b/Geo/MFace.h
index 9e7eb86730cc5921ba039b089744083b3c58d424..591351d6b43266c30e86442f7a98488b1790deb4 100644
--- a/Geo/MFace.h
+++ b/Geo/MFace.h
@@ -23,6 +23,9 @@ class MFace {
inline int getNumVertices() const { return _v[3] ? 4 : 3; }
inline MVertex *getVertex(const int i) const { return _v[i]; }
inline MVertex *getSortedVertex(const int i) const { return _v[int(_si[i])]; }
+
+ bool computeCorrespondence(const MFace&,int&,bool&) const;
+
void getOrderedVertices(std::vector<MVertex*> &verts) const
{
for(int i = 0; i < getNumVertices(); i++)
@@ -128,4 +131,7 @@ struct Less_Face : public std::binary_function<MFace, MFace, bool> {
}
};
+
+
+
#endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 638601d5242bd007c403e28780ba801239288705..5886e1f01e6b2d9623b1f8248b0226dee920a352 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -25,7 +25,8 @@ typedef std::map<std::pair<MVertex*,MVertex*>, std::vector<MVertex*> > edgeConta
// for each face (a list of vertices) we build a list of vertices that
// are the high order representation of the face
-typedef std::map<std::vector<MVertex*>, std::vector<MVertex*> > faceContainer;
+// typedef std::map<std::vector<MVertex*>, std::vector<MVertex*>> faceContainer;
+typedef std::map<MFace, std::vector<MVertex*>,Less_Face> faceContainer;
bool reparamOnFace(MVertex *v, GFace *gf, SPoint2 ¶m)
{
@@ -416,10 +417,16 @@ void getFaceVertices(GFace *gf,
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
- std::vector<MVertex*> p;
- face.getOrderedVertices(p);
- if(faceVertices.count(p)){
- vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
+ //std::vector<MVertex*> p;
+ // face.getOrderedVertices(p);
+ //if(faceVertices.count(p)){
+ // vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
+ // }
+
+ faceContainer::iterator fIter = faceVertices.find(face);
+
+ if(fIter!=faceVertices.end()){
+ vf.insert(vf.end(), fIter->second.begin(), fIter->second.end());
}
else{
SPoint2 pts[20];
@@ -432,6 +439,9 @@ void getFaceVertices(GFace *gf,
}
}
if(face.getNumVertices() == 3){ // triangles
+
+ std::vector<MVertex*>& vtcs = faceVertices[face];
+
for(int k = start ; k < points.size1() ; k++){
MVertex *v;
const double t1 = points(k, 0);
@@ -462,12 +472,17 @@ void getFaceVertices(GFace *gf,
v = new MVertex(X,Y,Z, gf);
}
}
- faceVertices[p].push_back(v);
+ // should be expensive -> induces a new search each time
+ // faceVertices[p].push_back(v);
+ vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
}
}
else if(face.getNumVertices() == 4){ // quadrangles
+
+ std::vector<MVertex*>& vtcs = faceVertices[face];
+
for(int j = 0; j < nPts; j++){
for(int k = 0; k < nPts; k++){
MVertex *v;
@@ -490,7 +505,8 @@ void getFaceVertices(GFace *gf,
GPoint pc = gf->point(uc, vc);
v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
}
- faceVertices[p].push_back(v);
+ // faceVertices[p].push_back(v);
+ vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
}
@@ -527,12 +543,21 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
- std::vector<MVertex*> p;
- face.getOrderedVertices(p);
- if(faceVertices.count(p)){
- vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
+ // std::vector<MVertex*> p;
+// face.getOrderedVertices(p);
+// if(faceVertices.count(p)){
+// vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
+// }
+
+ faceContainer::iterator fIter = faceVertices.find(face);
+
+ if (fIter!=faceVertices.end()) {
+ vf.insert(vf.end(), fIter->second.begin(), fIter->second.end());
}
else{
+
+ std::vector<MVertex*>& vtcs = faceVertices[face];
+
SPoint2 p0, p1, p2, p3;
bool reparamOK = true;
if(!linear &&
@@ -561,7 +586,8 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
v->setParameter (0,uc);
v->setParameter (1,vc);
}
- faceVertices[p].push_back(v);
+ // faceVertices[p].push_back(v);
+ vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
}
@@ -589,7 +615,8 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
GPoint pc = gf->point(uc, vc);
v = new MFaceVertex(pc.x(), pc.y(), pc.z(), gf, uc, vc);
}
- faceVertices[p].push_back(v);
+ // faceVertices[p].push_back(v);
+ vtcs.push_back(v);
gf->mesh_vertices.push_back(v);
vf.push_back(v);
}
@@ -599,25 +626,99 @@ void getFaceVertices(GFace *gf, MElement *ele, std::vector<MVertex*> &vf,
}
}
+void reorientTrianglePoints(std::vector<MVertex*>& vtcs,int orientation,bool swap) {
+
+
+ size_t nbPts = vtcs.size();
+
+ if (nbPts > 3) Msg::Error("Interior face nodes reorientation not supported for order > 4");
+ std::vector<MVertex*> tmp(nbPts);
+
+ // rotation
+ // --------
+
+ // --- interior "principal vertices"
+
+ // for (int i=0;i<3;i++) tmp[(i+orientation)%3] = vtcs[i];
+
+ for (int i=0;i<3;i++) tmp[i] = vtcs[(i+orientation)%3];
+
+
+ // normal swap
+ // -----------
+
+ if (swap) {
+ // --- interior "principal vertices"
+ for (int i=0;i<2;i++) for (int j=0;j<2-i;j++) vtcs[i*2+j] = tmp[i+j*2];
+ }
+
+ // no swap
+ // -------
+
+ else vtcs = tmp;
+}
+
+
+// KH: check face orientation wrt element ...
+
void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
faceContainer &faceVertices, bool linear, int nPts = 1)
{
+
for(int i = 0; i < ele->getNumFaces(); i++){
MFace face = ele->getFace(i);
- std::vector<MVertex*> p;
- face.getOrderedVertices(p);
- if(faceVertices.count(p)){
- vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
+ // std::vector<MVertex*> p;
+// face.getOrderedVertices(p);
+// if(faceVertices.count(p)){
+// // checking orientation not possible ???
+// vf.insert(vf.end(), faceVertices[p].begin(), faceVertices[p].end());
+// }
+
+ faceContainer::iterator fIter = faceVertices.find(face);
+
+ if (fIter!=faceVertices.end()) {
+
+ int orientation;
+ bool swap;
+
+ std::vector<MVertex*> vtcs = fIter->second;
+ if (fIter->first.computeCorrespondence(face,orientation,swap)) {
+ printf("element %d , face %d (%d,%d,%d) wrt (%d,%d,%d): orientation %d, swap %d\n",
+ ele->getNum(),i,
+ face.getVertex(0)->getNum(),
+ face.getVertex(1)->getNum(),
+ face.getVertex(2)->getNum(),
+ fIter->first.getVertex(0)->getNum(),
+ fIter->first.getVertex(1)->getNum(),
+ fIter->first.getVertex(2)->getNum(),
+ orientation, swap ? 1:0);
+ reorientTrianglePoints(vtcs,orientation,swap);
+ }
+ else Msg::Error("Error in face lookup for recuperation of high order face nodes");
+ vf.insert(vf.end(), vtcs.begin(), vtcs.end());
}
- else{
+ else{
+ std::vector<MVertex*>& vtcs = faceVertices[face];
+
if(face.getNumVertices() == 3){ // triangles
for(int j = 0; j < nPts; j++){
for(int k = 0 ; k < nPts - j - 1; k++){
- double t1 = (double)(j + 1) / (nPts + 1);
- double t2 = (double)(k + 1) / (nPts + 1);
+
+ // KH: inverted direction to stick with triangle vertex numbering
+ //
+ // 2
+ // | \
+ // 0 - 1
+ //
+ // double t1 = (double)(j + 1) / (nPts + 1);
+ // double t2 = (double)(k + 1) / (nPts + 1);
+ double t1 = (double)(k + 1) / (nPts + 1);
+ double t2 = (double)(j + 1) / (nPts + 1);
+
SPoint3 pc = face.interpolate(t1, t2);
MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
- faceVertices[p].push_back(v);
+ // faceVertices[p].push_back(v);
+ vtcs.push_back(v);
gr->mesh_vertices.push_back(v);
vf.push_back(v);
}
@@ -631,7 +732,8 @@ void getFaceVertices(GRegion *gr, MElement *ele, std::vector<MVertex*> &vf,
double t2 = 2. * (double)(k + 1) / (nPts + 1) - 1.;
SPoint3 pc = face.interpolate(t1, t2);
MVertex *v = new MVertex(pc.x(), pc.y(), pc.z(), gr);
- faceVertices[p].push_back(v);
+ // faceVertices[p].push_back(v);
+ vtcs.push_back(v);
gr->mesh_vertices.push_back(v);
vf.push_back(v);
}
diff --git a/Numeric/FunctionSpace.cpp b/Numeric/FunctionSpace.cpp
index 61b5597ec4111964955d248d0981f33e0f712fcc..40c6ce22d8008bfcee8bdff0e3aaa597c0a396a9 100644
--- a/Numeric/FunctionSpace.cpp
+++ b/Numeric/FunctionSpace.cpp
@@ -143,29 +143,34 @@ Double_Matrix generatePascalTetrahedron(int order)
int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
int nbdoftriangleserendip(int order) { return 3 * order; }
+//KH : caveat : node coordinates are not yet coherent with node numbering associated
+// to numbering of principal vertices of face !!!!
+
+// uv surface - orientation v0-v2-v1
void nodepositionface0(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= order; v[1]= 0; w[1] = 0.;
- u[2]= 0.; v[2]= order; w[2] = 0.;
+ u[1]= 0.; v[1]= order; w[1] = 0.;
+ u[2]= order; v[2]= 0.; w[2] = 0.;
// edges
for (int k = 0; k < (order - 1); k++){
- u[3 + k] = k + 1;
- v[3 + k] =0.;
+ u[3 + k] = 0.;
+ v[3 + k] = k + 1;
w[3 + k] = 0.;
- u[3 + order - 1 + k] = order - 1 - k ;
- v[3 + order - 1 + k] = k + 1;
+ u[3 + order - 1 + k] = k + 1;
+ v[3 + order - 1 + k] = order - 1 - k ;
w[3 + order - 1 + k] = 0.;
- u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ u[3 + 2 * (order - 1) + k] = order - 1 - k;
+ v[3 + 2 * (order - 1) + k] = 0.;
w[3 + 2 * (order - 1) + k] = 0.;
}
+
if (order > 2){
int nbdoftemp = nbdoftriangle(order - 3);
nodepositionface0(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
@@ -182,11 +187,11 @@ void nodepositionface0(int order, double *u, double *v, double *w)
w[k] = w[k] / order;
}
}
-
+// uw surface - orientation v0-v1-v3
void nodepositionface1(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
u[0]= 0.; v[0]= 0.; w[0] = 0.;
u[1]= order; v[1]= 0; w[1] = 0.;
@@ -222,34 +227,36 @@ void nodepositionface1(int order, double *u, double *v, double *w)
}
}
-void nodepositionface2(int order, double *u, double *v, double *w)
+// vw surface - orientation v0-v3-v2
+void nodepositionface2(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
- u[0]= order; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= 0.; v[2]= 0.; w[2] = order;
+ u[0]= 0.; v[0]= 0.; w[0] = 0.;
+ u[1]= 0.; v[1]= 0.; w[1] = order;
+ u[2]= 0.; v[2]= order; w[2] = 0.;
// edges
for (int k = 0; k < (order - 1); k++){
- u[3 + k] = order - 1 - k;
- v[3 + k] = k + 1;
- w[3 + k] = 0.;
+
+ u[3 + k] = 0.;
+ v[3 + k] = 0.;
+ w[3 + k] = k + 1;
u[3 + order - 1 + k] = 0.;
- v[3 + order - 1 + k] = order - 1 - k;
- w[3 + order - 1 + k] = k + 1;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = order - 1 - k;
- u[3 + 2 * (order - 1) + k] = k + 1;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = order - 1 - k;
+ u[3 + 2 * (order - 1) + k] = 0.;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = 0.;
}
if (order > 2){
int nbdoftemp = nbdoftriangle(order - 3);
nodepositionface2(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
&w[3 + 3 * (order - 1)]);
for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
}
@@ -261,34 +268,36 @@ void nodepositionface2(int order, double *u, double *v, double *w)
}
}
-void nodepositionface3(int order, double *u, double *v, double *w)
+// uvw surface - orientation v3-v1-v2
+void nodepositionface3(int order, double *u, double *v, double *w)
{
int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; return; }
+ if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
- u[0]= 0.; v[0]= 0.; w[0] = 0.;
- u[1]= 0.; v[1]= order; w[1] = 0.;
- u[2]= 0.; v[2]= 0.; w[2] = order;
+ u[0]= 0.; v[0]= 0.; w[0] = order;
+ u[1]= order; v[1]= 0.; w[1] = 0.;
+ u[2]= 0.; v[2]= order; w[2] = 0.;
// edges
for (int k = 0; k < (order - 1); k++){
- u[3 + k] = 0.;
- v[3 + k]= k + 1;
- w[3 + k] = 0.;
- u[3 + order - 1 + k] = 0.;
- v[3 + order - 1 + k] = order - 1 - k;
- w[3 + order - 1 + k] = k + 1;
+ u[3 + k] = k + 1;
+ v[3 + k] = 0.;
+ w[3 + k] = order - 1 - k;
+
+ u[3 + order - 1 + k] = order - 1 - k;
+ v[3 + order - 1 + k] = k + 1;
+ w[3 + order - 1 + k] = 0.;
u[3 + 2 * (order - 1) + k] = 0.;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = order - 1 - k;
+ v[3 + 2 * (order - 1) + k] = order - 1 - k;
+ w[3 + 2 * (order - 1) + k] = k + 1;
}
if (order > 2){
int nbdoftemp = nbdoftriangle(order - 3);
nodepositionface3(order - 3, &u[3 + 3 * (order - 1)], &v[3 + 3 * (order - 1)],
&w[3 + 3 * (order - 1)]);
for (int k = 0; k < nbdoftemp; k++){
- u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3);
+ u[3 + k + 3 * (order - 1)] = u[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
}
@@ -328,7 +337,8 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
point(3, 1) = 0.;
point(3, 2) = order;
- // edges e5 and e6 switched in original version
+ // edges e5 and e6 switched in original version, opposite direction
+ // the template has been defined in table edges_tetra and faces_tetra (MElement.h)
if (order > 1) {
for (int k = 0; k < (order - 1); k++) {
@@ -339,7 +349,7 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
// point(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
// point(4 + 5 * (order - 1) + k, 0) = 0.;
point(4 + 4 * (order - 1) + k, 0) = 0.;
- point(4 + 5 * (order - 1) + k, 0) = order - 1 - k;
+ point(4 + 5 * (order - 1) + k, 0) = k+1;
point(4 + k, 1) = 0.;
point(4 + order - 1 + k, 1) = k + 1;
@@ -347,15 +357,15 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
point(4 + 3 * (order - 1) + k, 1) = 0.;
// point(4 + 4 * (order - 1) + k, 1) = 0.;
// point(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
- point(4 + 4 * (order - 1) + k, 1) = order - 1 - k;
+ point(4 + 4 * (order - 1) + k, 1) = k+1;
point(4 + 5 * (order - 1) + k, 1) = 0.;
point(4 + k, 2) = 0.;
point(4 + order - 1 + k, 2) = 0.;
point(4 + 2 * (order - 1) + k, 2) = 0.;
- point(4 + 3 * (order - 1) + k, 2) = k + 1;
- point(4 + 4 * (order - 1) + k, 2) = k + 1;
- point(4 + 5 * (order - 1) + k, 2) = k + 1;
+ point(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
+ point(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
+ point(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
}
if (order > 2) {
@@ -368,6 +378,8 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
nodepositionface0(order - 3, u, v, w);
+ // u-v plane
+
for (int i = 0; i < nbdofface; i++){
point(ns + i, 0) = u[i] * (order - 3) + 1.;
point(ns + i, 1) = v[i] * (order - 3) + 1.;
@@ -376,6 +388,8 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
ns = ns + nbdofface;
+ // u-w plane
+
nodepositionface1(order - 3, u, v, w);
for (int i=0; i < nbdofface; i++){
@@ -383,23 +397,27 @@ Double_Matrix gmshGeneratePointsTetrahedron(int order, bool serendip)
point(ns + i, 1) = v[i] * (order - 3) ;
point(ns + i, 2) = w[i] * (order - 3) + 1.;
}
+
+ // v-w plane
ns = ns + nbdofface;
nodepositionface2(order - 3, u, v, w);
for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3) + 1.;
+ point(ns + i, 0) = u[i] * (order - 3);
point(ns + i, 1) = v[i] * (order - 3) + 1.;
point(ns + i, 2) = w[i] * (order - 3) + 1.;
}
+ // u-v-w plane
+
ns = ns + nbdofface;
nodepositionface3(order - 3, u, v, w);
for (int i = 0; i < nbdofface; i++){
- point(ns + i, 0) = u[i] * (order - 3);
+ point(ns + i, 0) = u[i] * (order - 3) + 1.;
point(ns + i, 1) = v[i] * (order - 3) + 1.;
point(ns + i, 2) = w[i] * (order - 3) + 1.;
}