diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 2900eaa8acda5f5a4270c976d16b5bea0533504f..2194c5a07c962bc2b1a14e25de0fe72b0d0e98a5 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1415,6 +1415,21 @@ MElement *MElementFactory::create(int type, std::vector<MVertex*> &v,
case MSH_PRI_6: return new MPrism(v, num, part);
case MSH_PRI_15: return new MPrism15(v, num, part);
case MSH_PRI_18: return new MPrism18(v, num, part);
+ case MSH_PRI_40: return new MPrismN(v, 3, num, part);
+ case MSH_PRI_75: return new MPrismN(v, 4, num, part);
+ case MSH_PRI_126: return new MPrismN(v, 5, num, part);
+ case MSH_PRI_196: return new MPrismN(v, 6, num, part);
+ case MSH_PRI_288: return new MPrismN(v, 7, num, part);
+ case MSH_PRI_405: return new MPrismN(v, 8, num, part);
+ case MSH_PRI_550: return new MPrismN(v, 9, num, part);
+ case MSH_PRI_24: return new MPrismN(v, 3, num, part);
+ case MSH_PRI_33: return new MPrismN(v, 4, num, part);
+ case MSH_PRI_42: return new MPrismN(v, 5, num, part);
+ case MSH_PRI_51: return new MPrismN(v, 6, num, part);
+ case MSH_PRI_60: return new MPrismN(v, 7, num, part);
+ case MSH_PRI_69: return new MPrismN(v, 8, num, part);
+ case MSH_PRI_78: return new MPrismN(v, 9, num, part);
+ case MSH_PRI_1: return new MPrismN(v, 0, num, part);
case MSH_TET_20: return new MTetrahedronN(v, 3, num, part);
case MSH_TET_35: return new MTetrahedronN(v, 4, num, part);
case MSH_TET_56: return new MTetrahedronN(v, 5, num, part);
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index beda50e53c21cf6f846fc010720235dbd20e8c48..3f11e25a5ea75b20db1fcab238572541227fe923 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -176,6 +176,10 @@ void MPrism18::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
_myGetEdgeRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
}
+void MPrismN::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;
}
@@ -184,6 +188,10 @@ int MPrism18::getNumEdgesRep() {
return 9 * CTX::instance()->mesh.numSubEdges;
}
+int MPrismN::getNumEdgesRep() {
+ return 9 * CTX::instance()->mesh.numSubEdges;
+}
+
void _myGetFaceRep(MPrism *pri, int num, double *x, double *y, double *z,
SVector3 *n, int numSubEdges)
@@ -413,6 +421,11 @@ void MPrism18::getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
_myGetFaceRep(this, num, x, y, z, n, CTX::instance()->mesh.numSubEdges);
}
+void MPrismN::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);
}
@@ -420,3 +433,88 @@ int MPrism15::getNumFacesRep() {
int MPrism18::getNumFacesRep() {
return 4 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
}
+
+int MPrismN::getNumFacesRep() {
+ return 4 * (CTX::instance()->mesh.numSubEdges * CTX::instance()->mesh.numSubEdges * 2);
+}
+
+namespace {
+
+void _addEdgeNodes(int num, bool reverse, int order, const std::vector<MVertex*> &vs,
+ int &ind, std::vector<MVertex*> &v)
+{
+
+ const int nNode = order-1, startNode = num*nNode, endNode = startNode+nNode-1;
+
+ if (reverse)
+ for (int i=endNode; i>=startNode; i--, ind++) v[ind] = vs[i];
+ else
+ for (int i=startNode; i<=endNode; i++, ind++) v[ind] = vs[i];
+
+}
+
+void _addFaceNodes(int num, int order, const std::vector<MVertex*> &vs,
+ int &ind, std::vector<MVertex*> &v)
+{
+
+ const int nNodeEd = order-1, nNodeTri = (order-2)*(order-1)/2;
+
+ int startNode, endNode;
+ if (num < 2) {
+ startNode = 9*nNodeEd+num*nNodeTri;
+ endNode = startNode+nNodeTri;
+ }
+ else {
+ const int nNodeQuad = (order-1)*(order-1);
+ startNode = 9*nNodeEd+2*nNodeTri+(num-2)*nNodeQuad;
+ endNode = startNode+nNodeQuad;
+ }
+
+ for (int i=startNode; i<endNode; i++, ind++) v[ind] = vs[i];
+
+}
+
+}
+
+// To be tested
+void MPrismN::getFaceVertices(const int num, std::vector<MVertex*> &v) const
+{
+
+ static const int edge[5][4] = {
+ {1, 3, 0, -1},
+ {6, 8, 7, -1},
+ {0, 4, 6, 2},
+ {2, 7, 5, 1},
+ {3, 5, 8, 4}
+ };
+ static const bool reverse[5][4] = {
+ {false, true, true, false},
+ {false, false, true, false},
+ {false, false, true, true},
+ {false, false, true, true},
+ {false, false, true, true}
+ };
+
+ int nNodeTotal, nEdge;
+ if (num < 2) { // Triangular face
+ nNodeTotal = (_order+1)*(_order+2)/2;
+ nEdge = 3;
+ }
+ else { // Quad face
+ nNodeTotal = (_order+1)*(_order+1);
+ nEdge = 4;
+ }
+
+ v.resize(nNodeTotal);
+
+ // Add corner nodes (there are nEdge corner nodes)
+ MPrism::_getFaceVertices(num, v);
+ int ind = nEdge;
+
+ // Add edge nodes
+ for (int iE=0; iE<nEdge; iE++) _addEdgeNodes(edge[num][iE],reverse[num][iE],_order,_vs,ind,v);
+
+ // Add face nodes
+ _addFaceNodes(num,_order,_vs,ind,v);
+
+}
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 964d8dd9d3056c37328ae200c057abed7825068a..2e27dc44979ac6065f67698857faaf05473d3e8e 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -385,4 +385,108 @@ class MPrism18 : public MPrism {
}
};
+/*
+ * MPrismN
+ */
+class MPrismN : public MPrism {
+ protected:
+ std::vector<MVertex *> _vs;
+ const char _order;
+ public :
+ MPrismN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3, MVertex *v4,
+ MVertex *v5, const std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MPrism(v0, v1, v2, v3, v4, v5, num, part), _vs(v), _order(order)
+ {
+ for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
+ }
+ MPrismN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
+ : MPrism(v, num, part), _order(order)
+ {
+ for (unsigned int i = 6; i < v.size(); i++) _vs.push_back(v[i]);
+ for (unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(2);
+ }
+ ~MPrismN(){}
+ virtual int getPolynomialOrder() const { return _order; }
+ virtual int getNumVertices() const { return 6+_vs.size(); }
+ virtual MVertex *getVertex(int num){ return num < 6 ? _v[num] : _vs[num-6]; }
+ virtual const MVertex *getVertex(int num) const{ return num < 6 ? _v[num] : _vs[num-6]; }
+ virtual int getNumEdgeVertices() const { return 9*(_order-1); }
+ virtual int getNumFaceVertices() const { int n = _order-1; return n*((n-1)+3*n); }
+ virtual int getNumVolumeVertices() const { int n = _order-1; return _vs.size()-n*(9+(n-1)+3*n); }
+ 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(_order+1);
+ MPrism::_getEdgeVertices(num, v);
+ const int n = _order-1;
+ for (int i=0; i<n; i++) v[2+i] = _vs[num*n+i];
+ }
+ 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;
+ virtual int getTypeForMSH() const {
+ switch (_order) {
+ case 0:
+ return MSH_PRI_1;
+ case 1:
+ return MSH_PRI_6;
+ case 2:
+ if (_vs.size() == 12) return MSH_PRI_18;
+ else if (_vs.size() == 9) return MSH_PRI_15;
+ break;
+ case 3:
+ if (_vs.size() == 34) return MSH_PRI_40;
+ else if (_vs.size() == 18) return MSH_PRI_24;
+ break;
+ case 4:
+ if (_vs.size() == 69) return MSH_PRI_75;
+ else if (_vs.size() == 27) return MSH_PRI_33;
+ break;
+ case 5:
+ if (_vs.size() == 120) return MSH_PRI_126;
+ else if (_vs.size() == 36) return MSH_PRI_42;
+ break;
+ case 6:
+ if (_vs.size() == 190) return MSH_PRI_196;
+ else if (_vs.size() == 45) return MSH_PRI_51;
+ break;
+ case 7:
+ if (_vs.size() == 282) return MSH_PRI_288;
+ else if (_vs.size() == 54) return MSH_PRI_60;
+ break;
+ case 8:
+ if (_vs.size() == 399) return MSH_PRI_405;
+ else if (_vs.size() == 63) return MSH_PRI_69;
+ break;
+ case 9:
+ if (_vs.size() == 544) return MSH_PRI_550;
+ else if (_vs.size() == 72) return MSH_PRI_78;
+ break;
+ }
+ Msg::Error("No tag matches a p%d prism with %d vertices", _order, 6+_vs.size());
+ return 0;
+ }
+ virtual const char *getStringForPOS() const {
+ switch (_order) {
+ case 0: return "SI0";
+ case 1: return "SI1";
+ case 2: return "SI2";
+ case 3: return "SI3";
+ case 4: return "SI4";
+ case 5: return "SI5";
+ case 6: return "SI6";
+ case 7: return "SI7";
+ case 8: return "SI8";
+ case 9: return "SI9";
+ }
+ return "";
+ }
+ virtual void getNode(int num, double &u, double &v, double &w) const
+ {
+ const fullMatrix<double> &p = getFunctionSpace()->points;
+ u = p(num,0); v = p(num,1); w = p(num,2);
+ }
+};
+
#endif
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index 69d6b8c09bd9379da24d79c0df670ecb0227d883..e8fd5ae08c6f12a608884ab5f4049f18f8a5de61 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -585,9 +585,7 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
switch (nPts){
case 0: return;
case 1: return;
- case 2:
- //BasisFactory::getNodalBasis(MSH_TET_20)->points.print();
- points = BasisFactory::getNodalBasis(MSH_TET_20)->points; break;
+ case 2: points = BasisFactory::getNodalBasis(MSH_TET_20)->points; break;
case 3: points = BasisFactory::getNodalBasis(MSH_TET_35)->points; break;
case 4: points = BasisFactory::getNodalBasis(MSH_TET_56)->points; break;
case 5: points = BasisFactory::getNodalBasis(MSH_TET_84)->points; break;
@@ -599,7 +597,7 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
}
- start = ((nPts+2)*(nPts+3)*(nPts+4)-(nPts-2)*(nPts-1)*(nPts))/6;
+ start = ((nPts+2)*(nPts+3)*(nPts+4)-(nPts-2)*(nPts-1)*(nPts))/6; // = 4+6*(p-1)+4*(p-2)*(p-1)/2 = 4*(p+1)*(p+2)/2-6*(p-1)-2*4 = 2*p*p+2
break;
case TYPE_HEX :
switch (nPts){
@@ -616,7 +614,24 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
}
- start = (nPts+2)*(nPts+2)*(nPts+2) - (nPts)*(nPts)*(nPts) ;
+ start = (nPts+2)*(nPts+2)*(nPts+2) - (nPts)*(nPts)*(nPts) ; // = 6*(p-1)*(p-1)+12*(p-1)+8 = 6*(p+1)*(p+1)-12*(p-1)-2*8 = 6*p*p+2
+ break;
+ case TYPE_PRI :
+ switch (nPts){
+ case 0: return;
+ case 1: return;
+ case 2: points = BasisFactory::getNodalBasis(MSH_PRI_40)->points; break;
+ case 3: points = BasisFactory::getNodalBasis(MSH_PRI_75)->points; break;
+ case 4: points = BasisFactory::getNodalBasis(MSH_PRI_126)->points; break;
+ case 5: points = BasisFactory::getNodalBasis(MSH_PRI_196)->points; break;
+ case 6: points = BasisFactory::getNodalBasis(MSH_PRI_288)->points; break;
+ case 7: points = BasisFactory::getNodalBasis(MSH_PRI_405)->points; break;
+ case 8: points = BasisFactory::getNodalBasis(MSH_PRI_550)->points; break;
+ default:
+ Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
+ break;
+ }
+ start = 4*(nPts+1)*(nPts+1)+2; // = 4*p*p+2 = 6+9*(p-1)+2*(p-2)*(p-1)/2+3*(p-1)*(p-1) = 2*(p+1)*(p+2)/2+3*(p+1)*(p+1)-9*(p-1)-2*6
break;
case TYPE_PYR:
switch (nPts){
@@ -857,8 +872,9 @@ static MPrism *setHighOrder(MPrism *p, GRegion *gr,
0, p->getPartition());
}
else{
- Msg::Error("PrismN generation not implemented");
- return 0;
+ return new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+ p->getVertex(3), p->getVertex(4), p->getVertex(5),
+ ve, nPts + 1, 0, p->getPartition());
}
}
else{
@@ -871,8 +887,20 @@ static MPrism *setHighOrder(MPrism *p, GRegion *gr,
0, p->getPartition());
}
else {
- Msg::Error("PrismN generation not implemented");
- return 0;
+ // create serendipity element to place internal vertices (we used to
+ // saturate face vertices also, but the corresponding function spaces do
+ // not exist anymore, and there is no theoretical justification for doing
+ // it either way)
+ MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+ p->getVertex(3), p->getVertex(4), p->getVertex(5),
+ ve, nPts + 1, 0, p->getPartition());
+ getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+ ve.insert(ve.end(), vf.begin(), vf.end());
+ getRegionVertices(gr, &incpl, p, vr, linear, nPts);
+ ve.insert(ve.end(), vr.begin(), vr.end());
+ return new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+ p->getVertex(3), p->getVertex(4), p->getVertex(5),
+ ve, nPts + 1, 0, p->getPartition());
}
}
}