diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index a20c54a3b87c52c40660fd82a0c5bb5bf2ddc833..70661ddc11a70881e4aba468b7ab8fde10c3e2df 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -244,7 +244,7 @@ int GmshBatch()
solver_batch_cb(0, (void*)CTX::instance()->launchSolverAtStartup);
#endif
- if (false) {
+ if (true) {
// 11/06/13 : This will be removed later !
for (int i = 1; i < MSH_NUM_TYPE+1; ++i) {
if (i == 76 || i == 77 || i == 78)
@@ -258,8 +258,10 @@ int GmshBatch()
continue;
}
if (MElement::ParentTypeFromTag(i) == TYPE_PYR) {
- Msg::Warning("ignoring pyramid %d", i);
- continue;
+ if (i > 124 && i < 132) {
+ Msg::Warning("ignoring serendipity pyramid %d (bug in Bergot implementation)", i);
+ continue;
+ }
}
if (MElement::ParentTypeFromTag(i) == TYPE_POLYG) {
Msg::Warning("ignoring polygone %d", i);
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 394720ec305f16ed6693d8e3e607e209c52d1def..2d7976a55cc851ef3d95890ad5cb239f8f692659 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1280,13 +1280,20 @@ int MElement::getInfoMSH(const int typeMSH, const char **const name)
case MSH_PYR_5 : if(name) *name = "Pyramid 5"; return 5;
case MSH_PYR_13 : if(name) *name = "Pyramid 13"; return 5 + 8;
case MSH_PYR_14 : if(name) *name = "Pyramid 14"; return 5 + 8 + 1;
- case MSH_PYR_30 : if(name) *name = "Pyramid 30"; return 5 + 8*2 + 4*1 + 1*4 + 1;
- case MSH_PYR_55 : if(name) *name = "Pyramid 55"; return 5 + 8*3 + 4*3 + 1*9 + 5;
- case MSH_PYR_91 : if(name) *name = "Pyramid 91"; return 91;
- case MSH_PYR_140 : if(name) *name = "Pyramid 140"; return 140;
- case MSH_PYR_204 : if(name) *name = "Pyramid 204"; return 204;
- case MSH_PYR_285 : if(name) *name = "Pyramid 285"; return 285;
- case MSH_PYR_385 : if(name) *name = "Pyramid 385"; return 385;
+ case MSH_PYR_30 : if(name) *name = "Pyramid 30"; return 5 + 8*2 + 4*1 + 1*4 + 1;
+ case MSH_PYR_55 : if(name) *name = "Pyramid 55"; return 5 + 8*3 + 4*3 + 1*9 + 5;
+ case MSH_PYR_91 : if(name) *name = "Pyramid 91"; return 5 + 8*4 + 4*6 + 1*16 + 14;
+ case MSH_PYR_140 : if(name) *name = "Pyramid 140"; return 5 + 8*5 + 4*10 + 1*25 + 30;
+ case MSH_PYR_204 : if(name) *name = "Pyramid 204"; return 5 + 8*6 + 4*15 + 1*36 + 55;
+ case MSH_PYR_285 : if(name) *name = "Pyramid 285"; return 5 + 8*7 + 4*21 + 1*49 + 91;
+ case MSH_PYR_385 : if(name) *name = "Pyramid 385"; return 5 + 8*8 + 4*28 + 1*64 + 140;
+ case MSH_PYR_29 : if(name) *name = "Pyramid 29"; return 5 + 8*2 + 4*1 + 1*4;
+ case MSH_PYR_50 : if(name) *name = "Pyramid 50"; return 5 + 8*3 + 4*3 + 1*9;
+ case MSH_PYR_77 : if(name) *name = "Pyramid 77"; return 5 + 8*4 + 4*6 + 1*16;
+ case MSH_PYR_110 : if(name) *name = "Pyramid 110"; return 5 + 8*5 + 4*10 + 1*25;
+ case MSH_PYR_149 : if(name) *name = "Pyramid 149"; return 5 + 8*6 + 4*15 + 1*36;
+ case MSH_PYR_194 : if(name) *name = "Pyramid 194"; return 5 + 8*7 + 4*21 + 1*49;
+ case MSH_PYR_245 : if(name) *name = "Pyramid 245"; return 5 + 8*8 + 4*28 + 1*64;
case MSH_POLYH_ : if(name) *name = "Polyhedron"; return 0;
case MSH_PNT_SUB : if(name) *name = "Point Xfem"; return 1;
case MSH_LIN_SUB : if(name) *name = "Line Xfem"; return 2;
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 465240b9088a225192ce1cbfd312cdc69c8ebc33..19fdbc628d28c72a00d4125d9b702472a889265e 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -166,11 +166,12 @@ class MPyramid : public MElement {
static int faces_pyramid(const int face, const int vert)
{
// only triangular faces
- static const int f[4][3] = {
- {0, 1, 4},
- {3, 0, 4},
- {1, 2, 4},
- {2, 3, 4}
+ static const int f[5][4] = {
+ {0, 1, 4, -1},
+ {3, 0, 4, -1},
+ {1, 2, 4, -1},
+ {2, 3, 4, -1},
+ {0, 3, 2, 1}
};
return f[face][vert];
}
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 21a69e2eb5dbcf7960ecbc44dfe8de5907dc7022..8fb9be5916ba7239fdb1e3ace21f95734f152546 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -1606,7 +1606,7 @@ nodalBasis::nodalBasis(int tag)
closureRef.resize(48, 0);
}
else {
- if (order < 2) generateFaceClosurePrism(closures, order);
+ generateFaceClosurePrism(closures, order);
generateFaceClosurePrismFull(fullClosures, closureRef, order);
}
break;
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 5ae836323f1400b65c0b7fd3808dd3d6ead67b16..4d12d058af3f39c868b4564399e5abd6ef369840 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -9,6 +9,7 @@
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
+#include "MPyramid.h"
/* --- Lines --- */
@@ -962,7 +963,8 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
{
- int nbMonomials = forSerendipPoints ? 6 + (order-1)*9 : (order + 1) * (order + 1)*(order + 2)/2;
+ int nbMonomials = forSerendipPoints ? 6 + (order-1)*9 :
+ (order + 1) * (order + 1)*(order + 2)/2;
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
@@ -1028,9 +1030,7 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
i2 = MPrism::faces_prism(iface, 2);
dudv.setAsProxy(dudvT);
}
- else {
- continue;
- }
+ else continue;
int u[3];
u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
@@ -1068,7 +1068,8 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints)
{
- int nbMonomials = forSerendipPoints ? 8 + (order-1)*12 : (order+1)*(order+1)*(order+1);
+ int nbMonomials = forSerendipPoints ? 8 + (order-1)*12 :
+ (order+1)*(order+1)*(order+1);
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
@@ -1121,10 +1122,10 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
}
}
- fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 2);
- dudv.add(1);
-
if (!forSerendipPoints) {
+ fullMatrix<int> dudv = gmshGenerateMonomialsQuadrangle(order - 2);
+ dudv.add(1);
+
for (int iface = 0; iface < 6; ++iface) {
int i0 = MHexahedron::faces_hexa(iface, 0);
int i1 = MHexahedron::faces_hexa(iface, 1);
@@ -1155,3 +1156,97 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
return monomials;
}
+fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints)
+{
+ int nbMonomials = forSerendipPoints ? 5 + (order-1)*8 :
+ (order+1)*((order+1)+1)*(2*(order+1)+1)/6;
+ fullMatrix<int> monomials(nbMonomials, 3);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+ monomials(0, 2) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = order;
+ monomials(1, 1) = 0;
+ monomials(1, 2) = 0;
+
+ monomials(2, 0) = order;
+ monomials(2, 1) = order;
+ monomials(2, 2) = 0;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = order;
+ monomials(3, 2) = 0;
+
+ monomials(4, 0) = 0;
+ monomials(4, 1) = 0;
+ monomials(4, 2) = order;
+
+ if (order > 1) {
+ int index = 5;
+ for (int iedge = 0; iedge < 8; ++iedge) {
+ int i0 = MPyramid::edges_pyramid(iedge, 0);
+ int i1 = MPyramid::edges_pyramid(iedge, 1);
+
+ int u_1 = (monomials(i1,0)-monomials(i0,0)) / order;
+ int u_2 = (monomials(i1,1)-monomials(i0,1)) / order;
+ int u_3 = (monomials(i1,2)-monomials(i0,2)) / order;
+
+ for (int i = 1; i < order; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + i * u_1;
+ monomials(index, 1) = monomials(i0, 1) + i * u_2;
+ monomials(index, 2) = monomials(i0, 2) + i * u_3;
+ }
+ }
+
+ if (!forSerendipPoints) {
+ fullMatrix<int> dudvQ = gmshGenerateMonomialsQuadrangle(order - 2);
+ dudvQ.add(1);
+
+ fullMatrix<int> dudvT;
+ if (order > 2) dudvT = gmshGenerateMonomialsTriangle(order - 3);
+ dudvT.add(1);
+
+ for (int iface = 0; iface < 5; ++iface) {
+ int i0, i1, i2;
+ i0 = MPyramid::faces_pyramid(iface, 0);
+ i1 = MPyramid::faces_pyramid(iface, 1);
+ fullMatrix<int> dudv;
+ if (MPyramid::faces_pyramid(iface, 3) != -1) {
+ i2 = MPyramid::faces_pyramid(iface, 3);
+ dudv.setAsProxy(dudvQ);
+ }
+ else if (order > 2) {
+ i2 = MPyramid::faces_pyramid(iface, 2);
+ dudv.setAsProxy(dudvT);
+ }
+ else continue;
+
+ int u[3];
+ u[0] = (monomials(i1, 0) - monomials(i0, 0)) / order;
+ u[1] = (monomials(i1, 1) - monomials(i0, 1)) / order;
+ u[2] = (monomials(i1, 2) - monomials(i0, 2)) / order;
+ int v[3];
+ v[0] = (monomials(i2, 0) - monomials(i0, 0)) / order;
+ v[1] = (monomials(i2, 1) - monomials(i0, 1)) / order;
+ v[2] = (monomials(i2, 2) - monomials(i0, 2)) / order;
+
+ for (int i = 0; i < dudv.size1(); ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + u[0] * dudv(i, 0) + v[0] * dudv(i, 1);
+ monomials(index, 1) = monomials(i0, 1) + u[1] * dudv(i, 0) + v[1] * dudv(i, 1);
+ monomials(index, 2) = monomials(i0, 2) + u[2] * dudv(i, 0) + v[2] * dudv(i, 1);
+ }
+ }
+
+ if (order > 2) {
+ fullMatrix<int> inner = gmshGenerateMonomialsPyramid(order - 3);
+ inner.add(1);
+ monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
+ }
+ }
+ }
+ }
+ return monomials;
+}
+
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index 4524266ac9288c2cf14d727818128dd20bd86397..9529cefff67313747e2686d3b7b8b1f53b648e73 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -44,7 +44,8 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
//fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order); //FIXME
-//fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool serendip);
+fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
+//fullMatrix<int> gmshGenerateMonomialsPyramidSerendipity(int order); //FIXME
diff --git a/Numeric/polynomialBasis.h b/Numeric/polynomialBasis.h
index c15e2eaed28a12507613e2505a256a7852f64c55..5c67c3580b3ef1f9548cbfed20be1bbefea358a1 100644
--- a/Numeric/polynomialBasis.h
+++ b/Numeric/polynomialBasis.h
@@ -75,8 +75,6 @@ class polynomialBasis : public nodalBasis
polynomialBasis(int tag);
~polynomialBasis();
- //int compareNewAlgoPointsWithOld() const;
-
virtual inline int getNumShapeFunctions() const {return coefficients.size1();}
virtual void f(double u, double v, double w, double *sf) const;
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index a977a06d6e5f0cd45d3cef115f7e2ec77e0a34a1..13c030b81b6c5d5c9e432d4f3fe855f9c05324e7 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -7,6 +7,15 @@
#include "pointsGenerators.h"
+static void copy(const fullMatrix<int> &orig, fullMatrix<double> &b)
+{
+ b.resize(orig.size1(), orig.size2());
+ for (int i = 0; i < orig.size1(); ++i) {
+ for (int j = 0; j < orig.size2(); ++j) {
+ b(i, j) = static_cast<double>(orig(i, j));
+ }
+ }
+}
pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
{
@@ -28,6 +37,21 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
delete[] fval;
+
+ // TEST NEW ALGO POINTS / MONOMIAL
+
+ monomials_newAlgo = gmshGenerateMonomialsPyramid(order, serendip);
+ copy(monomials_newAlgo, points_newAlgo);
+ if (order == 0) return;
+
+ for (int i = 0; i < points_newAlgo.size1(); ++i) {
+ points_newAlgo(i, 2) = points_newAlgo(i, 2) / order;
+ if (i != 4) {
+ const double duv = -1. + points_newAlgo(i, 2);
+ points_newAlgo(i, 0) = duv + points_newAlgo(i, 0) * 2. / order;
+ points_newAlgo(i, 1) = duv + points_newAlgo(i, 1) * 2. / order;
+ }
+ }
}
@@ -55,7 +79,7 @@ void pyramidalBasis::f(double u, double v, double w, double *val) const
val[i] = 0.;
for (int j = 0; j < N; j++) val[i] += VDMinv(i,j)*fval[j];
}
-
+
delete[] fval;
}
@@ -77,7 +101,7 @@ void pyramidalBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf)
for (int j = 0; j < N; j++) sf(iPt,i) += VDMinv(i,j)*fval[j];
}
}
-
+
delete[] fval;
}