diff --git a/Common/Gmsh.cpp b/Common/Gmsh.cpp
index c17b913e65caf5c6d276fdaf934bb1c5c4a642d1..4bb15f44193a3fbe09ede3e2c9f9289bb7182f11 100644
--- a/Common/Gmsh.cpp
+++ b/Common/Gmsh.cpp
@@ -246,40 +246,84 @@ int GmshBatch()
if (false) {
// 11/06/13 : This will be removed later !
+ std::vector<int> compare;
+ std::vector<int> test;
for (int i = 1; i < MSH_NUM_TYPE+1; ++i) {
if (i == 76 || i == 77 || i == 78)
continue;
+
+ const char **name = new const char*[1];
+ MElement::getInfoMSH(i, name);
+
if (i == 67 || i == 68 || i == 70) {
- Msg::Warning("ignoring unknown %d", i);
- continue;
- }
- if (MElement::ParentTypeFromTag(i) == TYPE_PRI && MElement::OrderFromTag(i) > 2) {
- Msg::Warning("ignoring prism %d (different implementation)", i);
+ Msg::Info("ignoring unknown %d", i);
continue;
}
- if (MElement::ParentTypeFromTag(i) == TYPE_PYR) {
- 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);
+ Msg::Info("ignoring polygone %d", i);
continue;
}
if (MElement::ParentTypeFromTag(i) == TYPE_POLYH) {
- Msg::Warning("ignoring polyhŹdre %d", i);
+ Msg::Info("ignoring polyhedre %d", i);
continue;
}
if (MElement::ParentTypeFromTag(i) == TYPE_XFEM) {
- Msg::Warning("ignoring xfem %d", i);
+ Msg::Info("ignoring xfem %d", i);
continue;
}
- const nodalBasis *fs = BasisFactory::getNodalBasis(i);
-
- if (fs && !fs->compareNewAlgoPointsWithOld() && !fs->serendip) {
- fs->compareNewAlgoBaseFunctionsWithOld();
+ if (MElement::ParentTypeFromTag(i) == TYPE_PRI && MElement::OrderFromTag(i) > 2) {
+ Msg::Warning("%d:ignoring %s (different node order)", i, *name);
+ test.push_back(i);
+ }
+ else if (MElement::ParentTypeFromTag(i) == TYPE_PYR && MElement::SerendipityFromTag(i) > 1) {
+ Msg::Warning("%d:ignoring Serendipity %s (bug in Bergot implementation & new Basis function not implemented)", i, *name);
+ }
+ else if (MElement::DimensionFromTag(i) == 3 && MElement::SerendipityFromTag(i) > 1) {
+ Msg::Warning("%d:ignoring 3D Serendipity : %s (different definition)", i, *name);
+ test.push_back(i);
+ }
+ else {
+ const nodalBasis *fs = BasisFactory::getNodalBasis(i);
+ if (fs && !fs->compareNewAlgoPointsWithOld()) {
+ compare.push_back(i);
+ }
+ }
+ }
+ Msg::Info(" ");
+ for (int parent = 1; parent < 12; ++parent) {
+ for (unsigned int i = 0; i < compare.size(); ++i) {
+ if (MElement::ParentTypeFromTag(compare[i]) == parent) {
+ if (parent == TYPE_PYR) {
+ const char **name = new const char*[1];
+ MElement::getInfoMSH(compare[i], name);
+ Msg::Warning("%d, %s: ignoring, old implementation much better and will be kept", compare[i], *name);
+ }
+ else {
+ const nodalBasis *fs = BasisFactory::getNodalBasis(compare[i]);
+ fs->compareNewAlgoBaseFunctionsWithOld();
+ }
+ }
+ }
+ }
+ Msg::Info(" ");
+ Msg::Info("---------------------------");
+ for (unsigned int i = 0; i < test.size(); ++i) {
+ const char **name = new const char*[1];
+ MElement::getInfoMSH(test[i], name);
+ //Msg::Info("%d, %s:", test[i], *name);
+ const nodalBasis *fs = BasisFactory::getNodalBasis(test[i]);
+ }
+ Msg::Info("---------------------------");
+ for (int parent = 1; parent < 12; ++parent) {
+ for (unsigned int i = 0; i < test.size(); ++i) {
+ if (MElement::ParentTypeFromTag(test[i]) == parent) {
+ const nodalBasis *fs = BasisFactory::getNodalBasis(test[i]);
+ /*const char **name = new const char*[1];
+ MElement::getInfoMSH(test[i], name);
+ Msg::Warning("%d, %s: untested !", test[i], *name);*/
+ fs->testNewAlgoBaseFunctions();
+ }
}
}
}
diff --git a/Numeric/nodalBasis.cpp b/Numeric/nodalBasis.cpp
index 5856ae2a099b5a7b357546aeac7c8534edd43ed1..2e1a2d0db543907ab0f20d3e95cefa428746cafe 100644
--- a/Numeric/nodalBasis.cpp
+++ b/Numeric/nodalBasis.cpp
@@ -15,7 +15,7 @@ int nodalBasis::compareNewAlgoPointsWithOld() const
const char **name = new const char*[1];
MElement::getInfoMSH(type, name);
if (points_newAlgo.size1() == 0) {
- Msg::Warning("%d: pas de points (%d, %d, %d) %s", type, parentType, serendip, order, *name);
+ Msg::Error("%d: pas de points (%d, %d, %d) %s", type, parentType, serendip, order, *name);
return 1;
}
if (points_newAlgo.size1() != points.size1()) {
@@ -92,23 +92,70 @@ int nodalBasis::compareNewAlgoBaseFunctionsWithOld() const
sumOne[1] += newVal[j];
}
sumError = std::sqrt(sumError/ndof);
- /*if (sumError > 1e-4) {
+ /*if (sumError > 1e-5) {
Msg::Error("(%.2f, %.2f, %.2f) -> fold=%.2e / fnew=%.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumOne[0], sumOne[1]);
//for (int j = 0; j < ndof; ++j) {
// Msg::Info("old %.3e vs %.3e new", oldVal[j], newVal[j]);
//}
}*/
- if (sumError > 1e-13) {
- if (type < 10) Msg::Warning(" f(%.2f, %.2f, %.2f) bad precision diff %.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
- else if (type < 100) Msg::Warning(" f(%.2f, %.2f, %.2f) bad precision diff %.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
- else Msg::Warning(" f(%.2f, %.2f, %.2f) bad precision diff %.2e", P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
+ if (sumError > 1e-5) {
+ Msg::Error("%d, %s: f(%.2f, %.2f, %.2f) bad precision diff %.2e", type, *name, P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
+ return 1;
+ }
+ else if (sumError > 1e-13) {
+ Msg::Warning("%d, %s: f(%.2f, %.2f, %.2f) bad precision diff %.2e", type, *name, P[0] / 1000., P[1] / 1000., P[2] / 1000., sumError);
+ return 1;
+ }
+ }
+ return 0;
+}
+
+int nodalBasis::testNewAlgoBaseFunctions() const
+{
+ const char **name = new const char*[1];
+ MElement::getInfoMSH(type, name);
+ int ndof = points_newAlgo.size1();
+ int P[3];
+ bool noproblem = true;
+ for (int i = 0; i < 10; ++i) {
+ if (i == 0) {
+ P[0] = 0;
+ P[1] = 0;
+ P[2] = 0;
+ }
+ else if (i == 1) {
+ P[0] = 0;
+ P[1] = 0;
+ P[2] = 1000;
+ }
+ else if (i == 2) {
+ P[0] = 1000;
+ P[1] = 0;
+ P[2] = 0;
+ }
+ else {
+ P[0] = std::rand() % 1001;
+ P[1] = std::rand() % (1001 - P[0]);
+ P[2] = std::rand() % (1001 - P[0] - P[1]);
+ }
+
+ double newVal[ndof];
+ fnew(P[0] / 1000., P[1] / 1000., P[2] / 1000., newVal);
+
+ double sumOne = 0;
+ for (int j = 0; j < ndof; ++j) {
+ sumOne += newVal[j];
+ }
+ if (sumOne > 1. + 1e-12 || sumOne < 1. - 1e-12) {
+ noproblem = false;
+ Msg::Warning("%d, %s: bad precision (%.2f, %.2f, %.2f) -> sum = %.2e (1 + %.2e)", type, *name, P[0] / 1000., P[1] / 1000., P[2] / 1000., sumOne, sumOne-1.);
+ //for (int j = 0; j < ndof; ++j) {
+ // Msg::Info(" %d : %.3e", j, newVal[j]);
+ //}
return 1;
}
- /*else if (type < 10) Msg::Info(" f(%.2f, %.2f, %.2f) ok", P[0] / 1000., P[1] / 1000., P[2] / 1000.);
- else if (type < 100) Msg::Info(" f(%.2f, %.2f, %.2f) ok", P[0] / 1000., P[1] / 1000., P[2] / 1000.);
- else Msg::Info(" f(%.2f, %.2f, %.2f) ok", P[0] / 1000., P[1] / 1000., P[2] / 1000.);
- */
}
+ if (noproblem) Msg::Info("%d, %s: no problem", type, *name);
return 0;
}
@@ -1477,7 +1524,10 @@ static void generateClosureOrder0(nodalBasis::clCont &closure, int nb)
nodalBasis::nodalBasis(int tag)
{
type = tag;
-
+ parentType = MElement::ParentTypeFromTag(tag);
+ order = MElement::OrderFromTag(tag);
+ serendip = MElement::SerendipityFromTag(tag) > 1;
+ /*
switch (tag) {
case MSH_PNT : parentType = TYPE_PNT; order = 0; serendip = false; break;
case MSH_LIN_1 : parentType = TYPE_LIN; order = 0; serendip = false; break;
@@ -1584,6 +1634,7 @@ nodalBasis::nodalBasis(int tag)
case MSH_HEX_296 : parentType = TYPE_HEX; order = 7; serendip = true; break;
case MSH_HEX_386 : parentType = TYPE_HEX; order = 8; serendip = true; break;
case MSH_HEX_488 : parentType = TYPE_HEX; order = 9; serendip = true; break;
+ case MSH_PYR_1 : parentType = TYPE_PYR; order = 0; serendip = false; break;
case MSH_PYR_5 : parentType = TYPE_PYR; order = 1; serendip = false; break;
case MSH_PYR_14 : parentType = TYPE_PYR; order = 2; serendip = false; break;
case MSH_PYR_30 : parentType = TYPE_PYR; order = 3; serendip = false; break;
@@ -1605,7 +1656,7 @@ nodalBasis::nodalBasis(int tag)
Msg::Error("Unknown function space %d: reverting to TET_4", tag);
parentType = TYPE_TET; order = 1; serendip = false; break;
}
-
+ */
switch (parentType) {
case TYPE_PNT :
diff --git a/Numeric/nodalBasis.h b/Numeric/nodalBasis.h
index c7339acca791daeb6f1f5bcdfeadc80bb98d0009..d9d48e911b79a7f8d0a7b5bf45d311cccb543da5 100644
--- a/Numeric/nodalBasis.h
+++ b/Numeric/nodalBasis.h
@@ -34,6 +34,7 @@ class nodalBasis {
int compareNewAlgoPointsWithOld() const;
int compareNewAlgoBaseFunctionsWithOld() const;
+ int testNewAlgoBaseFunctions() const;
// closures is the list of the nodes of each face, in the order of
// the polynomialBasis of the face; fullClosures is mapping of the
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index 20ff853be9a7ab268b727ef4d781de09f336fa0f..b56dea53126e5c53da73cf6943a30dfa8a3c353c 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -27,782 +27,103 @@ void getDoubleMonomials(fullMatrix<double>& doubleMon,
}
}
-/* --- Lines --- */
+// Points Generators
fullMatrix<double> gmshGeneratePointsLine(int order)
{
- fullMatrix<double> points(order + 1, 1);
- points(0,0) = 0;
- if (order > 0) {
- points(0, 0) = -1.;
- points(1, 0) = 1.;
- double dd = 2. / order;
- for (int i = 2; i < order + 1; i++) points(i, 0) = -1. + dd * (i - 1);
- }
+ fullMatrix<double> points;
+ getDoubleMonomials(points, gmshGenerateMonomialsLine, order, false);
+
+ if (order == 0) return points;
+
+ points.scale(2./order);
+ points.add(-1.);
return points;
}
-/* --- Triangles --- */
-
fullMatrix<double> gmshGeneratePointsTriangle(int order, bool serendip)
{
- int nbPoints = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
- fullMatrix<double> points(nbPoints, 2);
-
- points(0, 0) = 0;
- points(0, 1) = 0;
-
- if (order > 0) {
- double dd = 1. / order;
-
- points(1, 0) = 1;
- points(1, 1) = 0;
- points(2, 0) = 0;
- points(2, 1) = 1;
-
- int index = 3;
-
- if (order > 1) {
-
- double ksi = 0;
- double eta = 0;
-
- for (int i = 0; i < order - 1; i++, index++) {
- ksi += dd;
- points(index, 0) = ksi;
- points(index, 1) = eta;
- }
-
- ksi = 1.;
+ fullMatrix<double> points;
+ getDoubleMonomials(points, gmshGenerateMonomialsTriangle, order, serendip);
- for (int i = 0; i < order - 1; i++, index++) {
- ksi -= dd;
- eta += dd;
- points(index, 0) = ksi;
- points(index, 1) = eta;
- }
-
- eta = 1.;
- ksi = 0.;
-
- for (int i = 0; i < order - 1; i++, index++) {
- eta -= dd;
- points(index, 0) = ksi;
- points(index, 1) = eta;
- }
+ if (order == 0) return points;
- if (order > 2 && !serendip) {
- fullMatrix<double> inner = gmshGeneratePointsTriangle(order - 3, serendip);
- inner.scale(1. - 3. * dd);
- inner.add(dd);
- points.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
+ points.scale(1./order);
return points;
}
-/* --- Quadrangles --- */
-
fullMatrix<double> gmshGeneratePointsQuadrangle(int order, bool serendip)
{
- int nbPoints = serendip ? order*4 : (order+1)*(order+1);
- fullMatrix<double> points(nbPoints, 2);
+ fullMatrix<double> points;
+ getDoubleMonomials(points, gmshGenerateMonomialsQuadrangle, order, serendip);
- if (order > 0) {
- points(0, 0) = -1;
- points(0, 1) = -1;
- points(1, 0) = 1;
- points(1, 1) = -1;
- points(2, 0) = 1;
- points(2, 1) = 1;
- points(3, 0) = -1;
- points(3, 1) = 1;
+ if (order == 0) return points;
- if (order > 1) {
- int index = 4;
- const static int edges[4][2]={{0,1},{1,2},{2,3},{3,0}};
- for (int iedge=0; iedge<4; iedge++) {
- int p0 = edges[iedge][0];
- int p1 = edges[iedge][1];
- for (int i = 1; i < order; i++, index++) {
- points(index, 0) = points(p0, 0) + i*(points(p1,0)-points(p0,0))/order;
- points(index, 1) = points(p0, 1) + i*(points(p1,1)-points(p0,1))/order;
- }
- }
- if (order >= 2 && !serendip) {
- fullMatrix<double> inner = gmshGeneratePointsQuadrangle(order-2, false);
- inner.scale(1. - 2./order);
- points.copy(inner, 0, nbPoints - index, 0, 2, index, 0);
- }
- }
- }
- else {
- points(0, 0) = 0;
- points(0, 1) = 0;
- }
+ points.scale(2./order);
+ points.add(-1.);
return points;
-
-
-}
-
-/* --- Tetahedra --- */
-
-static int nbdoftriangle(int order) { return (order + 1) * (order + 2) / 2; }
-
-//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
-static void nodepositionface0(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- 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]= order; v[2]= 0.; 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] = k + 1;
- v[3 + order - 1 + k] = order - 1 - k ;
- w[3 + order - 1 + k] = 0.;
-
- 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)],
- &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.;
- 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);
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
}
-// uw surface - orientation v0-v1-v3
-static void nodepositionface1(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- 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]= 0.; w[2] = order;
- // edges
- for (int k = 0; k < (order - 1); k++){
- u[3 + k] = k + 1;
- v[3 + k] = 0.;
- w[3 + k] = 0.;
-
- u[3 + order - 1 + k] = order - 1 - k;
- v[3 + order - 1 + k] = 0.;
- w[3 + order - 1+ k ] = k + 1;
-
- u[3 + 2 * (order - 1) + k] = 0. ;
- v[3 + 2 * (order - 1) + k] = 0.;
- w[3 + 2 * (order - 1) + k] = order - 1 - k;
- }
- if (order > 2){
- int nbdoftemp = nbdoftriangle(order - 3);
- nodepositionface1(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.;
- v[3 + k + 3 * (order - 1)] = v[3 + k + 3 * (order - 1)] * (order - 3);
- w[3 + k + 3 * (order - 1)] = w[3 + k + 3 * (order - 1)] * (order - 3) + 1.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-// vw surface - orientation v0-v3-v2
-static 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]= 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] = 0.;
- v[3 + k] = 0.;
- w[3 + k] = k + 1;
-
- u[3 + order - 1 + k] = 0.;
- v[3 + order - 1 + k] = k + 1;
- w[3 + 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);
- 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.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-// uvw surface - orientation v3-v1-v2
-static void nodepositionface3(int order, double *u, double *v, double *w)
-{
- int ndofT = nbdoftriangle(order);
- if (order == 0) { u[0] = 0.; v[0] = 0.; w[0] = 0.; return; }
-
- 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] = 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] = 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) + 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.;
- }
- }
- for (int k = 0; k < ndofT; k++){
- u[k] = u[k] / order;
- v[k] = v[k] / order;
- w[k] = w[k] / order;
- }
-}
-
-
fullMatrix<double> gmshGeneratePointsTetrahedron(int order, bool serendip)
{
+ fullMatrix<double> points;
+ getDoubleMonomials(points, gmshGenerateMonomialsTetrahedron, order, serendip);
- int nbPoints =
- (serendip ?
- 4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
- (order + 1) * (order + 2) * (order + 3) / 6);
-
- fullMatrix<double> points(nbPoints, 3);
-
- double overOrder = (order == 0 ? 1. : 1. / order);
-
- points(0, 0) = 0.;
- points(0, 1) = 0.;
- points(0, 2) = 0.;
+ if (order == 0) return points;
- if (order > 0) {
- points(1, 0) = order;
- points(1, 1) = 0;
- points(1, 2) = 0;
-
- points(2, 0) = 0.;
- points(2, 1) = order;
- points(2, 2) = 0.;
-
- points(3, 0) = 0.;
- points(3, 1) = 0.;
- points(3, 2) = order;
-
- // 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++) {
- points(4 + k, 0) = k + 1;
- points(4 + order - 1 + k, 0) = order - 1 - k;
- points(4 + 2 * (order - 1) + k, 0) = 0.;
- points(4 + 3 * (order - 1) + k, 0) = 0.;
- // points(4 + 4 * (order - 1) + k, 0) = order - 1 - k;
- // points(4 + 5 * (order - 1) + k, 0) = 0.;
- points(4 + 4 * (order - 1) + k, 0) = 0.;
- points(4 + 5 * (order - 1) + k, 0) = k+1;
-
- points(4 + k, 1) = 0.;
- points(4 + order - 1 + k, 1) = k + 1;
- points(4 + 2 * (order - 1) + k, 1) = order - 1 - k;
- points(4 + 3 * (order - 1) + k, 1) = 0.;
- // points(4 + 4 * (order - 1) + k, 1) = 0.;
- // points(4 + 5 * (order - 1) + k, 1) = order - 1 - k;
- points(4 + 4 * (order - 1) + k, 1) = k+1;
- points(4 + 5 * (order - 1) + k, 1) = 0.;
-
- points(4 + k, 2) = 0.;
- points(4 + order - 1 + k, 2) = 0.;
- points(4 + 2 * (order - 1) + k, 2) = 0.;
- points(4 + 3 * (order - 1) + k, 2) = order - 1 - k;
- points(4 + 4 * (order - 1) + k, 2) = order - 1 - k;
- points(4 + 5 * (order - 1) + k, 2) = order - 1 - k;
- }
-
- if (order > 2) {
- int ns = 4 + 6 * (order - 1);
- int nbdofface = nbdoftriangle(order - 3);
-
- double *u = new double[nbdofface];
- double *v = new double[nbdofface];
- double *w = new double[nbdofface];
-
- nodepositionface0(order - 3, u, v, w);
-
- // u-v plane
-
- for (int i = 0; i < nbdofface; i++){
- points(ns + i, 0) = u[i] * (order - 3) + 1.;
- points(ns + i, 1) = v[i] * (order - 3) + 1.;
- points(ns + i, 2) = w[i] * (order - 3);
- }
-
- ns = ns + nbdofface;
-
- // u-w plane
-
- nodepositionface1(order - 3, u, v, w);
-
- for (int i=0; i < nbdofface; i++){
- points(ns + i, 0) = u[i] * (order - 3) + 1.;
- points(ns + i, 1) = v[i] * (order - 3) ;
- points(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++){
- points(ns + i, 0) = u[i] * (order - 3);
- points(ns + i, 1) = v[i] * (order - 3) + 1.;
- points(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++){
- points(ns + i, 0) = u[i] * (order - 3) + 1.;
- points(ns + i, 1) = v[i] * (order - 3) + 1.;
- points(ns + i, 2) = w[i] * (order - 3) + 1.;
- }
-
- ns = ns + nbdofface;
-
- delete [] u;
- delete [] v;
- delete [] w;
-
- if (!serendip && order > 3) {
-
- fullMatrix<double> interior = gmshGeneratePointsTetrahedron(order - 4, false);
- for (int k = 0; k < interior.size1(); k++) {
- points(ns + k, 0) = 1. + interior(k, 0) * (order - 4);
- points(ns + k, 1) = 1. + interior(k, 1) * (order - 4);
- points(ns + k, 2) = 1. + interior(k, 2) * (order - 4);
- }
- }
- }
- }
- }
-
- points.scale(overOrder);
+ points.scale(1./order);
return points;
-
}
-
-
-/* --- Hexahedra --- */
-
fullMatrix<double> gmshGeneratePointsHexahedron(int order, bool serendip)
{
- int nbPoints = (order+1)*(order+1)*(order+1);
- if (serendip) nbPoints -= (order-1)*(order-1)*(order-1);
- fullMatrix<double> points(nbPoints, 3);
-
- // should be a public member of MHexahedron, just copied
- static const int edges[12][2] = {
- {0, 1},
- {0, 3},
- {0, 4},
- {1, 2},
- {1, 5},
- {2, 3},
- {2, 6},
- {3, 7},
- {4, 5},
- {4, 7},
- {5, 6},
- {6, 7}
- };
- static const int faces[6][4] = {
- {0, 3, 2, 1},
- {0, 1, 5, 4},
- {0, 4, 7, 3},
- {1, 2, 6, 5},
- {2, 3, 7, 6},
- {4, 5, 6, 7}
- };
+ fullMatrix<double> points;
+ getDoubleMonomials(points, gmshGenerateMonomialsHexahedron, order, serendip);
- if (order > 0) {
- points(0, 0) = -1;
- points(0, 1) = -1;
- points(0, 2) = -1;
- points(1, 0) = 1;
- points(1, 1) = -1;
- points(1, 2) = -1;
- points(2, 0) = 1;
- points(2, 1) = 1;
- points(2, 2) = -1;
- points(3, 0) = -1;
- points(3, 1) = 1;
- points(3, 2) = -1;
-
- points(4, 0) = -1;
- points(4, 1) = -1;
- points(4, 2) = 1;
- points(5, 0) = 1;
- points(5, 1) = -1;
- points(5, 2) = 1;
- points(6, 0) = 1;
- points(6, 1) = 1;
- points(6, 2) = 1;
- points(7, 0) = -1;
- points(7, 1) = 1;
- points(7, 2) = 1;
+ if (order == 0) return points;
- if (order > 1) {
- int index = 8;
- for (int iedge=0; iedge<12; iedge++) {
- int p0 = edges[iedge][0];
- int p1 = edges[iedge][1];
- for (int i = 1; i < order; i++, index++) {
- points(index, 0) = points(p0, 0) + i*(points(p1,0)-points(p0,0))/order;
- points(index, 1) = points(p0, 1) + i*(points(p1,1)-points(p0,1))/order;
- points(index, 2) = points(p0, 2) + i*(points(p1,2)-points(p0,2))/order;
- }
- }
-
- fullMatrix<double> fp = gmshGeneratePointsQuadrangle(order - 2, false);
- fp.scale(1. - 2./order);
- for (int iface=0; iface<6; iface++) {
- int p0 = faces[iface][0];
- int p1 = faces[iface][1];
- int p2 = faces[iface][2];
- int p3 = faces[iface][3];
- for (int i=0;i<fp.size1();i++, index++){
- const double u = fp(i,0);
- const double v = fp(i,1);
- const double newU =
- 0.25*(1.-u)*(1.-v)*points(p0,0) +
- 0.25*(1.+u)*(1.-v)*points(p1,0) +
- 0.25*(1.+u)*(1.+v)*points(p2,0) +
- 0.25*(1.-u)*(1.+v)*points(p3,0) ;
- const double newV =
- 0.25*(1.-u)*(1.-v)*points(p0,1) +
- 0.25*(1.+u)*(1.-v)*points(p1,1) +
- 0.25*(1.+u)*(1.+v)*points(p2,1) +
- 0.25*(1.-u)*(1.+v)*points(p3,1) ;
- const double newW =
- 0.25*(1.-u)*(1.-v)*points(p0,2) +
- 0.25*(1.+u)*(1.-v)*points(p1,2) +
- 0.25*(1.+u)*(1.+v)*points(p2,2) +
- 0.25*(1.-u)*(1.+v)*points(p3,2) ;
- points(index, 0) = newU;
- points(index, 1) = newV;
- points(index, 2) = newW;
- }
- }
- if (!serendip) {
- fullMatrix<double> inner = gmshGeneratePointsHexahedron(order - 2, false);
- inner.scale(1. - 2./order);
- points.copy(inner, 0, nbPoints - index, 0, 3, index, 0);
- }
- }
- }
- else if (order == 0){
- points(0, 0) = 0;
- points(0, 1) = 0;
- points(0, 2) = 0;
- }
+ points.scale(2./order);
+ points.add(-1.);
return points;
}
-/* --- Prisms --- */
-
fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip)
{
- const double prism18Pts[18][3] = {
- {0, 0, -1}, // 0
- {1, 0, -1}, // 1
- {0, 1, -1}, // 2
- {0, 0, 1}, // 3
- {1, 0, 1}, // 4
- {0, 1, 1}, // 5
- {0.5, 0, -1}, // 6
- {0, 0.5, -1}, // 7
- {0, 0, 0}, // 8
- {0.5, 0.5, -1}, // 9
- {1, 0, 0}, // 10
- {0, 1, 0}, // 11
- {0.5, 0, 1}, // 12
- {0, 0.5, 1}, // 13
- {0.5, 0.5, 1}, // 14
- {0.5, 0, 0}, // 15
- {0, 0.5, 0}, // 16
- {0.5, 0.5, 0}, // 17
- };
-
- int nbPoints = (order + 1)*(order + 1)*(order + 2)/2;
- fullMatrix<double> points(nbPoints, 3);
-
- int index = 0;
- fullMatrix<double> triPoints = gmshGeneratePointsTriangle(order,false);
- fullMatrix<double> linePoints = gmshGeneratePointsLine(order);
-
- if (order == 2)
- for (int i =0; i<18; i++)
- for (int j=0; j<3;j++)
- points(i,j) = prism18Pts[i][j];
- else
- for (int j = 0; j <linePoints.size1() ; j++) {
- for (int i = 0; i < triPoints.size1(); i++) {
- points(index,0) = triPoints(i,0);
- points(index,1) = triPoints(i,1);
- points(index,2) = linePoints(j,0);
- index ++;
- }
- }
+ fullMatrix<double> points;
+ getDoubleMonomials(points, gmshGenerateMonomialsPrism, order, serendip);
- return points;
+ if (order == 0) return points;
-}
+ fullMatrix<double> tmp;
+ tmp.setAsProxy(points, 0, 2);
+ tmp.scale(1./order);
-/* --- Pyramids --- */
+ tmp.setAsProxy(points, 2, 1);
+ tmp.scale(2./order);
+ tmp.add(-1.);
+
+ return points;
+}
fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
{
- int nbPoints = (order+1)*((order+1)+1)*(2*(order+1)+1)/6;
-
- // Remove volume points if incomplete basis
- if (serendip && order > 2) nbPoints -= (order-2)*((order-2)+1)*(2*(order-2)+1)/6;
-
- fullMatrix<double> points(nbPoints, 3);
-
- static const int edges[8][2] = {
- {0, 1},
- {0, 3},
- {0, 4},
- {1, 2},
- {1, 4},
- {2, 3},
- {2, 4},
- {3, 4},
- };
- static const int faces[5][4] = {
- {0, 1, 4, -1},
- {3, 0, 4, -1},
- {1, 2, 4, -1},
- {2, 3, 4, -1},
- {0, 3, 2, 1}
- };
-
- if (order == 0) {
- points(0,0) = 0.0;
- points(0,1) = 0.0;
- points(0,2) = 0.0;
- return points;
- }
-
- // Principal vertices of the pyramid
- points(0,0) = -1.0; points(0,1) = -1.0; points(0,2) = 0.0;
- points(1,0) = 1.0; points(1,1) = -1.0; points(1,2) = 0.0;
- points(2,0) = 1.0; points(2,1) = 1.0; points(2,2) = 0.0;
- points(3,0) = -1.0; points(3,1) = 1.0; points(3,2) = 0.0;
- points(4,0) = 0.0; points(4,1) = 0.0; points(4,2) = 1.0;
-
- if (order > 1) {
- int p = 5;
-
- // Edges
- for (int e = 0; e < 8; e++) {
- double vec[3];
- vec[0] = points(edges[e][1],0) - points(edges[e][0],0);
- vec[1] = points(edges[e][1],1) - points(edges[e][0],1);
- vec[2] = points(edges[e][1],2) - points(edges[e][0],2);
- for (int i = 0; i < order - 1; i++) {
- points(p,0) = points(edges[e][0],0) + vec[0]/order * (i+1);
- points(p,1) = points(edges[e][0],1) + vec[1]/order * (i+1);
- points(p,2) = points(edges[e][0],2) + vec[2]/order * (i+1);
- p++;
- }
- }
-
- // Faces
- for (int f = 0; f < 4; f++) {
- fullMatrix<double> fp = gmshGeneratePointsTriangle(order, false);
-
- fullVector<double> u(4);
- fullVector<double> v(4);
- fullVector<double> w(4);
- fullVector<double> y(4);
-
- fullVector<double> up(4);
- fullVector<double> vp(4);
- fullVector<double> wp(4);
- fullVector<double> yp(4);
-
- for (int c = 0; c < 3; c++) {
- u(c) = fp(0,c);
- v(c) = fp(1,c);
- w(c) = fp(2,c);
- up(c) = points(faces[f][0],c);
- vp(c) = points(faces[f][1],c);
- wp(c) = points(faces[f][2],c);
- }
-
- u(2) = 0.0;
- v(2) = 0.0;
- w(2) = 0.0;
- u(3) = 1.0;
- v(3) = 1.0;
- w(3) = 1.0;
-
- y(0) = u(0) + ((v(1)-u(1))*(w(2)-u(2)) - (v(2)-u(2))*(w(1)-u(1)));
- y(1) = u(1) + ((v(2)-u(2))*(w(0)-u(0)) - (v(0)-u(0))*(w(2)-u(2)));
- y(2) = u(2) + ((v(0)-u(0))*(w(1)-u(1)) - (v(1)-u(1))*(w(0)-u(0)));
- y(3) = 1.0;
-
- up(3) = 1.0;
- vp(3) = 1.0;
- wp(3) = 1.0;
-
- yp(0) = up(0) + ((vp(1)-up(1))*(wp(2)-up(2)) - (vp(2)-up(2))*(wp(1)-up(1)));
- yp(1) = up(1) + ((vp(2)-up(2))*(wp(0)-up(0)) - (vp(0)-up(0))*(wp(2)-up(2)));
- yp(2) = up(2) + ((vp(0)-up(0))*(wp(1)-up(1)) - (vp(1)-up(1))*(wp(0)-up(0)));
- yp(3) = 1.0;
-
- fullMatrix<double> M(4,4);
- fullMatrix<double> Mp(4,4);
+ fullMatrix<double> points;
+ getDoubleMonomials(points, gmshGenerateMonomialsPyramid, order, serendip);
- M(0,0) = u(0); M(0,1) = v(0); M(0,2) = w(0); M(0,3) = y(0);
- M(1,0) = u(1); M(1,1) = v(1); M(1,2) = w(1); M(1,3) = y(1);
- M(2,0) = u(2); M(2,1) = v(2); M(2,2) = w(2); M(2,3) = y(2);
- M(3,0) = u(3); M(3,1) = v(3); M(3,2) = w(3); M(3,3) = y(3);
+ if (order == 0) return points;
- Mp(0,0) = up(0); Mp(0,1) = vp(0); Mp(0,2) = wp(0); Mp(0,3) = yp(0);
- Mp(1,0) = up(1); Mp(1,1) = vp(1); Mp(1,2) = wp(1); Mp(1,3) = yp(1);
- Mp(2,0) = up(2); Mp(2,1) = vp(2); Mp(2,2) = wp(2); Mp(2,3) = yp(2);
- Mp(3,0) = up(3); Mp(3,1) = vp(3); Mp(3,2) = wp(3); Mp(3,3) = yp(3);
-
-
- M.invertInPlace();
- fullMatrix<double> T(4,4);
- Mp.mult(M, T);
-
- for(int k = 3*order; k < fp.size1(); k++) {
- fullVector<double> pnt(4);
- pnt(0) = fp(k,0);
- pnt(1) = fp(k,1);
- pnt(2) = 0.0;
- pnt(3) = 1.0;
-
- fullVector<double> res(4);
- T.mult(pnt, res);
-
- points(p, 0) = res(0);
- points(p, 1) = res(1);
- points(p, 2) = res(2);
-
- p++;
-
- }
- }
-
- fullMatrix<double> fpq = gmshGeneratePointsQuadrangle(order-2, false);
- fullMatrix<double> transform(2,2);
- fullMatrix<double> scaled(fpq.size1(), fpq.size2());
- transform(0,0) = (order-2)/(double)order; transform(0,1) = 0.0;
- transform(1,0) = 0.0; transform(1,1) = (order-2)/(double)order;
- fpq.mult(transform, scaled);
-
- for (int i = 0; i < scaled.size1(); i++) {
- points(p, 0) = scaled(i, 1);
- points(p, 1) = scaled(i, 0);
- points(p, 2) = 0.0;
- p++;
- }
-
- // Volume
- if (!serendip && order > 2) {
- fullMatrix<double> volume_points = gmshGeneratePointsPyramid(order - 3, false);
- // scale to order-3/order
- fullMatrix<double> T(3,3);
- T(0,0) = (order - 3.0) / order; T(0,1) = 0.0; T(0,2) = 0.0;
- T(1,0) = 0.0; T(1,1) = (order - 3.0) / order; T(1,2) = 0.0;
- T(2,0) = 0.0; T(2,1) = 0.0; T(2,2) = (order - 3.0) / order;
- fullMatrix<double> scaled(volume_points.size1(), volume_points.size2());
- volume_points.mult(T, scaled);
-
- // translate 1/order
- for (int i = 0; i < scaled.size1(); i++) {
- points(p, 0) = scaled(i,0);
- points(p, 1) = scaled(i,1);
- points(p, 2) = scaled(i,2) + 1.0/order;
- p++;
- }
- }
- }
- return points;
+ for (int i = 0; i < points.size1(); ++i) {
+ points(i, 2) = points(i, 2) / order;
+ const double duv = -1. + points(i, 2);
+ points(i, 0) = duv + points(i, 0) * 2. / order;
+ points(i, 1) = duv + points(i, 1) * 2. / order;
+ }
+ return points;
}
-fullMatrix<int> gmshGenerateMonomialsLine(int order)
+// Monomials Generators
+
+fullMatrix<int> gmshGenerateMonomialsLine(int order, bool serendip)
{
fullMatrix<int> monomials(order + 1, 1);
monomials(0,0) = 0;
@@ -816,6 +137,7 @@ fullMatrix<int> gmshGenerateMonomialsLine(int order)
fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
{
int nbMonomials = serendip ? 3 * order : (order + 1) * (order + 2) / 2;
+ if (serendip && !order) nbMonomials = 1;
fullMatrix<int> monomials(nbMonomials, 2);
monomials(0, 0) = 0;
@@ -855,6 +177,7 @@ fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip)
fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? order*4 : (order+1)*(order+1);
+ if (forSerendipPoints && !order) nbMonomials = 1;
fullMatrix<int> monomials(nbMonomials, 2);
monomials(0, 0) = 0;
@@ -895,6 +218,44 @@ fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoint
return monomials;
}
+fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order)
+{
+ int nbMonomials = order ? order*4 : 1;
+ fullMatrix<int> monomials(nbMonomials, 2);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = 1;
+ monomials(1, 1) = 0;
+
+ monomials(2, 0) = 1;
+ monomials(2, 1) = 1;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = 1;
+
+ if (order > 1) {
+ int index = 3;
+ for (int p = 2; p <= order; ++p) {
+ monomials(++index, 0) = p;
+ monomials( index, 1) = 0;
+
+ monomials(++index, 0) = p;
+ monomials( index, 1) = 1;
+
+ monomials(++index, 0) = 1;
+ monomials( index, 1) = p;
+
+ monomials(++index, 0) = 0;
+ monomials( index, 1) = p;
+ }
+ }
+ }
+ return monomials;
+}
+
//KH : caveat : node coordinates are not yet coherent with node numbering associated
// to numbering of principal vertices of face !!!!
fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
@@ -904,6 +265,7 @@ fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip)
4 + 6 * std::max(0, order - 1) + 4 * std::max(0, (order - 2) * (order - 1) / 2) :
(order + 1) * (order + 2) * (order + 3) / 6);
+ if (serendip && !order) nbMonomials = 1;
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
@@ -983,6 +345,7 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 6 + (order-1)*9 :
(order + 1) * (order + 1)*(order + 2)/2;
+ if (forSerendipPoints && !order) nbMonomials = 1;
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
@@ -1084,10 +447,84 @@ fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints)
return monomials;
}
+fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order)
+{
+ int nbMonomials = order ? 6 + (order-1) * 9 : 1;
+ fullMatrix<int> monomials(nbMonomials, 3);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+ monomials(0, 2) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = 1;
+ monomials(1, 1) = 0;
+ monomials(1, 2) = 0;
+
+ monomials(2, 0) = 0;
+ monomials(2, 1) = 1;
+ monomials(2, 2) = 0;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = 0;
+ monomials(3, 2) = 1;
+
+ monomials(4, 0) = 1;
+ monomials(4, 1) = 0;
+ monomials(4, 2) = 1;
+
+ monomials(5, 0) = 0;
+ monomials(5, 1) = 1;
+ monomials(5, 2) = 1;
+
+ if (order > 1) {
+ const int ind[7][3] = {
+ {2, 0, 0},
+ {2, 0, 1},
+
+ {0, 2, 0},
+ {0, 2, 1},
+
+ {0, 0, 2},
+ {1, 0, 2},
+ {0, 1, 2}
+ };
+ int val[3] = {0, 1, -1};
+ int index = 5;
+ for (int p = 2; p <= order; ++p) {
+ val[2] = p;
+ for (int i = 0; i < 7; ++i) {
+ monomials(++index, 0) = val[ind[i][0]];
+ monomials( index, 1) = val[ind[i][1]];
+ monomials( index, 2) = val[ind[i][2]];
+ }
+ }
+
+ int val0 = 1;
+ int val1 = order - 1;
+ for (int p = 2; p <= order; ++p) {
+ monomials(++index, 0) = val0;
+ monomials( index, 1) = val1;
+ monomials( index, 2) = 0;
+
+ monomials(++index, 0) = val0;
+ monomials( index, 1) = val1;
+ monomials( index, 2) = 1;
+
+ ++val0;
+ --val1;
+ }
+ }
+ }
+ return monomials;
+}
+
+
fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints)
{
int nbMonomials = forSerendipPoints ? 8 + (order-1)*12 :
(order+1)*(order+1)*(order+1);
+ if (forSerendipPoints && !order) nbMonomials = 1;
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
@@ -1174,10 +611,83 @@ fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoint
return monomials;
}
+
+fullMatrix<int> gmshGenerateMonomialsHexaSerendipity(int order)
+{
+ int nbMonomials = order ? 8 + (order-1) * 12 : 1;
+ fullMatrix<int> monomials(nbMonomials, 3);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+ monomials(0, 2) = 0;
+
+ if (order > 0) {
+ monomials(1, 0) = 1;
+ monomials(1, 1) = 0;
+ monomials(1, 2) = 0;
+
+ monomials(2, 0) = 1;
+ monomials(2, 1) = 1;
+ monomials(2, 2) = 0;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = 1;
+ monomials(3, 2) = 0;
+
+ monomials(4, 0) = 0;
+ monomials(4, 1) = 0;
+ monomials(4, 2) = 1;
+
+ monomials(5, 0) = 1;
+ monomials(5, 1) = 0;
+ monomials(5, 2) = 1;
+
+ monomials(6, 0) = 1;
+ monomials(6, 1) = 1;
+ monomials(6, 2) = 1;
+
+ monomials(7, 0) = 0;
+ monomials(7, 1) = 1;
+ monomials(7, 2) = 1;
+
+ if (order > 1) {
+ const int ind[12][3] = {
+ {2, 0, 0},
+ {2, 0, 1},
+ {2, 1, 1},
+ {2, 1, 0},
+
+ {0, 2, 0},
+ {0, 2, 1},
+ {1, 2, 1},
+ {1, 2, 0},
+
+ {0, 0, 2},
+ {0, 1, 2},
+ {1, 1, 2},
+ {1, 0, 2}
+ };
+ int val[3] = {0, 1, -1};
+ int index = 7;
+ for (int p = 2; p <= order; ++p) {
+ val[2] = p;
+ for (int i = 0; i < 12; ++i) {
+ monomials(++index, 0) = val[ind[i][0]];
+ monomials( index, 1) = val[ind[i][1]];
+ monomials( index, 2) = val[ind[i][2]];
+ }
+ }
+ }
+ }
+ 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;
+ if (forSerendipPoints && !order) nbMonomials = 1;
fullMatrix<int> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index 9529cefff67313747e2686d3b7b8b1f53b648e73..dc47ba71897d108fa54e2961a7ef974c48e2e0d4 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -32,17 +32,17 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip);
// Monomial exponents
-fullMatrix<int> gmshGenerateMonomialsLine(int order);
+fullMatrix<int> gmshGenerateMonomialsLine(int order, bool serendip = false);
fullMatrix<int> gmshGenerateMonomialsTriangle(int order, bool serendip = false);
fullMatrix<int> gmshGenerateMonomialsQuadrangle(int order, bool forSerendipPoints = false);
-//fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order); //FIXME
+fullMatrix<int> gmshGenerateMonomialsQuadSerendipity(int order);
fullMatrix<int> gmshGenerateMonomialsTetrahedron(int order, bool serendip = false);
fullMatrix<int> gmshGenerateMonomialsHexahedron(int order, bool forSerendipPoints = false);
-//fullMatrix<int> gmshGenerateMonomialsHexaSerendipity(int order); //FIXME
+fullMatrix<int> gmshGenerateMonomialsHexaSerendipity(int order);
fullMatrix<int> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints = false);
-//fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order); //FIXME
+fullMatrix<int> gmshGenerateMonomialsPrismSerendipity(int order);
fullMatrix<int> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
//fullMatrix<int> gmshGenerateMonomialsPyramidSerendipity(int order); //FIXME
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 50dfccdf826bbc320942e982a4510b846871a1af..5a7b626b9b19ced3fbdd8f3121678f92b4a081c5 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -433,7 +433,7 @@ static fullMatrix<double> generateLagrangeMonomialCoefficients
max = std::max(max, std::abs(unity(i, j)));
}
}
- if (max > 1e-10) Msg::Info(" max unity = %.3e", max);
+ //if (max > 1e-10) Msg::Info(" max unity = %.3e", max);
return coefficient;
}
@@ -481,60 +481,47 @@ polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
coefficients = generateLagrangeMonomialCoefficients(monomials, points);
- // TEST NEW ALGO POINTS / MONOMIALS
+ // NEW ALGO POINTS / MONOMIALS
- int rescale = 0;
switch (parentType) {
case TYPE_PNT :
+ points_newAlgo = gmshGeneratePointsLine(0);
monomials_newAlgo = gmshGenerateMonomialsLine(0);
break;
case TYPE_LIN :
+ points_newAlgo = gmshGeneratePointsLine(order);
monomials_newAlgo = gmshGenerateMonomialsLine(order);
- rescale = 1;
break;
case TYPE_TRI :
+ points_newAlgo = gmshGeneratePointsTriangle(order, serendip);
monomials_newAlgo = gmshGenerateMonomialsTriangle(order, serendip);
break;
case TYPE_QUA :
- monomials_newAlgo = gmshGenerateMonomialsQuadrangle(order, serendip);
- rescale = 1;
+ points_newAlgo = gmshGeneratePointsQuadrangle(order, serendip);
+ if (!serendip)
+ monomials_newAlgo = gmshGenerateMonomialsQuadrangle(order);
+ else
+ monomials_newAlgo = gmshGenerateMonomialsQuadSerendipity(order);
break;
case TYPE_TET :
+ points_newAlgo = gmshGeneratePointsTetrahedron(order, serendip);
monomials_newAlgo = gmshGenerateMonomialsTetrahedron(order, serendip);
break;
case TYPE_HEX :
- monomials_newAlgo = gmshGenerateMonomialsHexahedron(order, serendip);
- rescale = 1;
+ points_newAlgo = gmshGeneratePointsHexahedron(order, serendip);
+ if (!serendip)
+ monomials_newAlgo = gmshGenerateMonomialsHexahedron(order);
+ else
+ monomials_newAlgo = gmshGenerateMonomialsHexaSerendipity(order);
break;
case TYPE_PRI :
- monomials_newAlgo = gmshGenerateMonomialsPrism(order, serendip);
- rescale = 2;
+ points_newAlgo = gmshGeneratePointsPrism(order, serendip);
+ if (!serendip)
+ monomials_newAlgo = gmshGenerateMonomialsPrism(order);
+ else
+ monomials_newAlgo = gmshGenerateMonomialsPrismSerendipity(order);
break;
}
- copy(monomials_newAlgo, points_newAlgo);
- //if (order == 0) return;
- switch (rescale) {
- case 0 :
- points_newAlgo.scale(1./order);
- break;
- case 1 :
- points_newAlgo.scale(2./order);
- points_newAlgo.add(-1.);
- break;
- case 2 :
- {
- fullMatrix<double> tmp;
- tmp.setAsProxy(points_newAlgo, 0, 2);
- tmp.scale(1./order);
-
- tmp.setAsProxy(points_newAlgo, 2, 1);
- tmp.scale(2./order);
- tmp.add(-1.);
- break;
- }
- default :
- break;
- }
fullMatrix<double> monDouble;
switch (parentType) {
@@ -551,6 +538,7 @@ polynomialBasis::polynomialBasis(int tag) : nodalBasis(tag)
/*else*/ copy(monomials_newAlgo, monDouble);
break;
}
+
coefficients_newAlgo = generateLagrangeMonomialCoefficients(monDouble, points_newAlgo);
}
diff --git a/Numeric/pyramidalBasis.cpp b/Numeric/pyramidalBasis.cpp
index 888752d749e907a55c55fd94f7d54bd4cf7ddc2d..1ee2fd4a0a0b9111e89c4091a0c94e53e17cbeac 100644
--- a/Numeric/pyramidalBasis.cpp
+++ b/Numeric/pyramidalBasis.cpp
@@ -62,7 +62,7 @@ static fullMatrix<double> generateLagrangeMonomialCoefficientsPyr
max = std::max(max, std::abs(unity(i, j)));
}
}
- if (max > 1e-10) Msg::Info("mon max unity = %.3e", max);
+ //if (max > 1e-14) Msg::Info("mon max unity = %.3e", max);
return coefficient;
}
@@ -98,23 +98,13 @@ pyramidalBasis::pyramidalBasis(int tag) : nodalBasis(tag)
max = std::max(max, std::abs(unity(i, j)));
}
}
- if (max > 1e-10) Msg::Info("leg max unity = %.3e", max);
+ //if (max > 1e-14) Msg::Info("leg max unity = %.3e", max);
// 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;
- }
- }
+ points_newAlgo = gmshGeneratePointsPyramid(order, serendip);
fullMatrix<double> monDouble;
copy(monomials_newAlgo, monDouble);