diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 18a16f4536f844a9530416305ad924b37126a197..f364c55674a05e1c43b27b5a237485a80b72c53c 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -6,6 +6,7 @@
#include "MPrism.h"
#include "Numeric.h"
#include "BasisFactory.h"
+#include "Context.h"
int MPrism::getVolumeSign()
{
@@ -159,3 +160,295 @@ void MPrism::getFaceInfo(const MFace &face, int &ithFace, int &sign, int &rot) c
}
Msg::Error("Could not get face information for prism %d", getNum());
}
+
+static void _myGetEdgeRep(MPrism *pri, int num, double *x, double *y, double *z,
+ SVector3 *n, int numSubEdges) {
+
+ //const int numSubEdges = CTX::instance()->mesh.numSubEdges;
+ static double pp[6][3] = {
+ {0,0,-1},{1,0,-1},{0,1,-1},
+ {0,0,1},{1,0,1},{0,1,1} };
+ static int ed [9][2] = {
+ {0,1},{0,2},{0,3},{1,2},{1,4},{2,5},
+ {3,4},{3,5},{4,5}
+ };
+ int iEdge = num / numSubEdges;
+ int iSubEdge = num % numSubEdges;
+
+ int iVertex1 = ed [iEdge][0];
+ int iVertex2 = ed [iEdge][1];
+ double t1 = (double) iSubEdge / (double) numSubEdges;
+ double u1 = pp[iVertex1][0] * (1.-t1) + pp[iVertex2][0] * t1;
+ double v1 = pp[iVertex1][1] * (1.-t1) + pp[iVertex2][1] * t1;
+ double w1 = pp[iVertex1][2] * (1.-t1) + pp[iVertex2][2] * t1;
+
+ double t2 = (double) (iSubEdge+1) / (double) numSubEdges;
+ double u2 = pp[iVertex1][0] * (1.-t2) + pp[iVertex2][0] * t2;
+ double v2 = pp[iVertex1][1] * (1.-t2) + pp[iVertex2][1] * t2;
+ double w2 = pp[iVertex1][2] * (1.-t2) + pp[iVertex2][2] * t2;
+
+ SPoint3 pnt1, pnt2;
+ pri->pnt(u1, v1, w1, pnt1);
+ pri->pnt(u2, v2, w2, pnt2);
+ x[0] = pnt1.x(); x[1] = pnt2.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z();
+
+ // not great, but better than nothing
+ //static const int f[6] = {0, 0, 0, 1, 2, 3};
+ n[0] = n[1] = 1 ;
+}
+
+
+void MPrism15::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n) {
+ _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+void MPrism18::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n) {
+ _myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+int MPrism15::getNumEdgesRep() {
+ return 9 * CTX::instance()->mesh.numSubEdges;
+}
+
+int MPrism18::getNumEdgesRep() {
+ return 9 * CTX::instance()->mesh.numSubEdges;
+}
+
+
+void _myGetFaceRep(MPrism *pri, int num, double *x, double *y, double *z,
+ SVector3 *n, int numSubEdges)
+{
+ static double pp[6][3] = {
+ {0,0,-1},{1,0,-1},{0,1,-1},
+ {0,0,1},{1,0,1},{0,1,1} };
+
+ int iFace = num / (numSubEdges * numSubEdges);
+ int iSubFace = num % (numSubEdges * numSubEdges);
+
+ if (iFace > 1) {
+ iFace = num / (2*numSubEdges * numSubEdges) + 1;
+ iSubFace = num % (2*numSubEdges * numSubEdges);
+ }
+
+ int iVertex1 = pri->faces_prism(iFace,0);
+ int iVertex2 = pri->faces_prism(iFace,1);
+ int iVertex3 = pri->faces_prism(iFace,2);
+ int iVertex4 = pri->faces_prism(iFace,3);
+
+ SPoint3 pnt1, pnt2, pnt3;
+ // double J1[3][3], J2[3][3], J3[3][3];
+
+ /*
+ 0
+ 0 1
+ 0 1 2
+ 0 1 2 3
+ 0 1 2 3 4
+ 0 1 2 3 4 5
+ */
+
+ // on each layer, we have (numSubEdges) * 2 triangles
+ // ix and iy are the coordinates of the sub-quadrangle
+
+ if (iFace > 1) {
+ int io = iSubFace%2;
+ int ix = (iSubFace/2)/numSubEdges;
+ int iy = (iSubFace/2)%numSubEdges;
+
+ const double d = 2. / numSubEdges;
+ double ox = -1. + d*ix;
+ double oy = -1. + d*iy;
+
+ if (io == 0){
+ double U1 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V1 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W1 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ ox += d;
+
+ double U2 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V2 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W2 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ oy += d;
+
+ double U3 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V3 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W3 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ pri->pnt(U1, V1, W1, pnt1);
+ pri->pnt(U2, V2, W2, pnt2);
+ pri->pnt(U3, V3, W3, pnt3);
+ }
+ else{
+ double U1 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V1 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W1 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ ox += d;
+ oy += d;
+
+ double U2 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V2 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W2 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ ox -= d;
+
+ double U3 =
+ pp[iVertex1][0] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][0] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][0] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][0] * (1.-ox)*(1+oy)*.25;
+ double V3 =
+ pp[iVertex1][1] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][1] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][1] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][1] * (1.-ox)*(1+oy)*.25;
+ double W3 =
+ pp[iVertex1][2] * (1.-ox)*(1-oy)*.25 +
+ pp[iVertex2][2] * (1.+ox)*(1-oy)*.25 +
+ pp[iVertex3][2] * (1.+ox)*(1+oy)*.25 +
+ pp[iVertex4][2] * (1.-ox)*(1+oy)*.25;
+
+ pri->pnt(U1, V1, W1, pnt1);
+ pri->pnt(U2, V2, W2, pnt2);
+ pri->pnt(U3, V3, W3, pnt3);
+ }
+
+ } else {
+ int ix = 0, iy = 0;
+ int nbt = 0;
+ for (int i = 0; i < numSubEdges; i++){
+ int nbl = (numSubEdges - i - 1) * 2 + 1;
+ nbt += nbl;
+ if (nbt > iSubFace){
+ iy = i;
+ ix = nbl - (nbt - iSubFace);
+ break;
+ }
+ }
+
+ const double d = 1. / numSubEdges;
+
+ double u1, v1, u2, v2, u3, v3;
+ if (ix % 2 == 0){
+ u1 = ix / 2 * d; v1= iy*d;
+ u2 = (ix / 2 + 1) * d ; v2 = iy * d;
+ u3 = ix / 2 * d ; v3 = (iy+1) * d;
+ }
+ else{
+ u1 = (ix / 2 + 1) * d; v1= iy * d;
+ u2 = (ix / 2 + 1) * d; v2= (iy + 1) * d;
+ u3 = ix / 2 * d ; v3 = (iy + 1) * d;
+ }
+
+ double U1 = pp[iVertex1][0] * (1.-u1-v1) + pp[iVertex2][0] * u1 + pp[iVertex3][0] * v1;
+ double U2 = pp[iVertex1][0] * (1.-u2-v2) + pp[iVertex2][0] * u2 + pp[iVertex3][0] * v2;
+ double U3 = pp[iVertex1][0] * (1.-u3-v3) + pp[iVertex2][0] * u3 + pp[iVertex3][0] * v3;
+
+ double V1 = pp[iVertex1][1] * (1.-u1-v1) + pp[iVertex2][1] * u1 + pp[iVertex3][1] * v1;
+ double V2 = pp[iVertex1][1] * (1.-u2-v2) + pp[iVertex2][1] * u2 + pp[iVertex3][1] * v2;
+ double V3 = pp[iVertex1][1] * (1.-u3-v3) + pp[iVertex2][1] * u3 + pp[iVertex3][1] * v3;
+
+ double W1 = pp[iVertex1][2] * (1.-u1-v1) + pp[iVertex2][2] * u1 + pp[iVertex3][2] * v1;
+ double W2 = pp[iVertex1][2] * (1.-u2-v2) + pp[iVertex2][2] * u2 + pp[iVertex3][2] * v2;
+ double W3 = pp[iVertex1][2] * (1.-u3-v3) + pp[iVertex2][2] * u3 + pp[iVertex3][2] * v3;
+
+ pri->pnt(U1, V1, W1, pnt1);
+ pri->pnt(U2, V2, W2, pnt2);
+ pri->pnt(U3, V3, W3, pnt3);
+ }
+
+
+
+ x[0] = pnt1.x(); x[1] = pnt2.x(); x[2] = pnt3.x();
+ y[0] = pnt1.y(); y[1] = pnt2.y(); y[2] = pnt3.y();
+ z[0] = pnt1.z(); z[1] = pnt2.z(); z[2] = pnt3.z();
+
+ SVector3 d1(x[1] - x[0], y[1] - y[0], z[1] - z[0]);
+ SVector3 d2(x[2] - x[0], y[2] - y[0], z[2] - z[0]);
+ n[0] = crossprod(d1, d2);
+ n[0].normalize();
+ n[1] = n[0];
+ n[2] = n[0];
+}
+
+void MPrism15::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+void MPrism18::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
+{
+ _myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
+}
+
+int MPrism15::getNumFacesRep() {
+ return 4 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
+}
+
+int MPrism18::getNumFacesRep() {
+ return 4 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
+}
\ No newline at end of file
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 618e10fb30cbd704196dd90e0dbdc9e6a9b654b7..964d8dd9d3056c37328ae200c057abed7825068a 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -241,42 +241,16 @@ class MPrism15 : public MPrism {
}
virtual MVertex *getVertexINP(int num){ return getVertexBDF(num); }
virtual int getNumEdgeVertices() const { return 9; }
- virtual int getNumEdgesRep(){ return 18; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[18][2] = {
- {0, 6}, {6, 1},
- {0, 7}, {7, 2},
- {0, 8}, {8, 3},
- {1, 9}, {9, 2},
- {1, 10}, {10, 4},
- {2, 11}, {11, 5},
- {3, 12}, {12, 4},
- {3, 13}, {13, 5},
- {4, 14}, {14, 5}
- };
- static const int f[9] = {0, 1, 2, 0, 2, 3, 1, 1, 1};
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
- }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MPrism::_getEdgeVertices(num, v);
v[2] = _vs[num];
}
- virtual int getNumFacesRep(){ return 26; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[26][3] = {
- {0, 7, 6}, {2, 9, 7}, {1, 6, 9}, {6, 7, 9},
- {3, 12, 13}, {4, 14, 12}, {5, 13, 14}, {12, 14, 13},
- {0, 6, 8}, {1, 10, 6}, {4, 12, 10}, {3, 8, 12}, {6, 10, 12}, {6, 12, 8},
- {0, 8, 7}, {3, 13, 8}, {5, 11, 13}, {2, 7, 11}, {7, 8, 13}, {7, 13, 11},
- {1, 9, 10}, {2, 11, 9}, {5, 14, 11}, {4, 10, 14}, {9, 11, 14}, {9, 14, 10}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
+ virtual int getNumFacesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize((num < 2) ? 6 : 8);
@@ -361,45 +335,16 @@ class MPrism18 : public MPrism {
virtual const MVertex *getVertex(int num) const{ return num < 6 ? _v[num] : _vs[num - 6]; }
virtual int getNumEdgeVertices() const { return 9; }
virtual int getNumFaceVertices() const { return 3; }
- virtual int getNumEdgesRep(){ return 18; }
- virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int e[18][2] = {
- {0, 6}, {6, 1},
- {0, 7}, {7, 2},
- {0, 8}, {8, 3},
- {1, 9}, {9, 2},
- {1, 10}, {10, 4},
- {2, 11}, {11, 5},
- {3, 12}, {12, 4},
- {3, 13}, {13, 5},
- {4, 14}, {14, 5}
- };
- static const int f[9] = {0, 1, 2, 0, 2, 3, 1, 1, 1};
- _getEdgeRep(getVertex(e[num][0]), getVertex(e[num][1]), x, y, z, n, f[num / 2]);
- }
+ virtual int getNumEdgesRep();
+ virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getEdgeVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize(3);
MPrism::_getEdgeVertices(num, v);
v[2] = _vs[num];
}
- virtual int getNumFacesRep(){ return 32; }
- virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
- {
- static const int f[32][3] = {
- {0, 7, 6}, {2, 9, 7}, {1, 6, 9}, {6, 7, 9},
- {3, 12, 13}, {4, 14, 12}, {5, 13, 14}, {12, 14, 13},
- {0, 6, 15}, {0, 15, 8}, {1, 10, 15}, {1, 15, 6},
- {4, 12, 15}, {4, 15, 10}, {3, 8, 15}, {3, 15, 12},
- {0, 8, 16}, {0, 16, 7}, {3, 13, 16}, {3, 16, 8},
- {5, 11, 16}, {5, 16, 13}, {2, 7, 16}, {2, 16, 11},
- {1, 9, 17}, {1, 17, 10}, {2, 11, 17}, {2, 17, 9},
- {5, 14, 17}, {5, 17, 11}, {4, 10, 17}, {4, 17, 14}
- };
- _getFaceRep(getVertex(f[num][0]), getVertex(f[num][1]), getVertex(f[num][2]),
- x, y, z, n);
- }
+ virtual int getNumFacesRep();
+ virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n);
virtual void getFaceVertices(const int num, std::vector<MVertex*> &v) const
{
v.resize((num < 2) ? 6 : 9);