diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index e23bb8ca36893c3e27357e0a2cd2001a2c7d24ff..402cddd299ffd7bd3f5a361d2dd0fcf31f40d3fe 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -676,7 +676,7 @@ void JacobianBasis::interpolate(const fullVector<double> &jacobian,
fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
{
- int nbMonomials = (order+3)*((order+3)+1)*(2*(order+3)+1)/6 - 5;
+ const int nbMonomials = (order+3)*(order+3)*(order+1);
fullMatrix<double> monomials(nbMonomials, 3);
if (order == 0) {
@@ -706,28 +706,28 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
monomials(4, 1) = 0;
monomials(4, 2) = order;
- monomials(5, 0) = 2;
+ monomials(5, 0) = order+2;
monomials(5, 1) = 0;
monomials(5, 2) = order;
- monomials(6, 0) = 2;
- monomials(6, 1) = 2;
+ monomials(6, 0) = order+2;
+ monomials(6, 1) = order+2;
monomials(6, 2) = order;
monomials(7, 0) = 0;
- monomials(7, 1) = 2;
+ monomials(7, 1) = order+2;
monomials(7, 2) = order;
int index = 8;
- static const int bottom_edges[8][2] = {
+ static const int bottom_edges[4][2] = {
{0, 1},
{1, 2},
{2, 3},
{3, 0}
};
- // bottom "edges"
+ // bottom & top "edges"
for (int iedge = 0; iedge < 4; ++iedge) {
int i0 = bottom_edges[iedge][0];
int i1 = bottom_edges[iedge][1];
@@ -740,22 +740,14 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
monomials(index, 1) = monomials(i0, 1) + i * u_2;
monomials(index, 2) = 0;
}
+ for (int i = 1; i < order + 2; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + i * u_1;
+ monomials(index, 1) = monomials(i0, 1) + i * u_2;
+ monomials(index, 2) = order;
+ }
}
- // top "edges"
- for (int iedge = 0; iedge < 4; ++iedge, ++index) {
- int i0 = bottom_edges[iedge][0] + 4;
- int i1 = bottom_edges[iedge][1] + 4;
-
- int u_1 = (monomials(i1,0)-monomials(i0,0)) / 2;
- int u_2 = (monomials(i1,1)-monomials(i0,1)) / 2;
-
- monomials(index, 0) = monomials(i0, 0) + u_1;
- monomials(index, 1) = monomials(i0, 1) + u_2;
- monomials(index, 2) = monomials(i0, 2);
- }
-
- // bottom "face"
+ // bottom & top "face"
fullMatrix<double> uv = gmshGenerateMonomialsQuadrangle(order);
uv.add(1);
for (int i = 0; i < uv.size1(); ++i, ++index) {
@@ -763,16 +755,15 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
monomials(index, 1) = uv(i, 1);
monomials(index, 2) = 0;
}
-
- // top "face"
- monomials(index, 0) = 1;
- monomials(index, 1) = 1;
- monomials(index, 2) = order;
- ++index;
+ for (int i = 0; i < uv.size1(); ++i, ++index) {
+ monomials(index, 0) = uv(i, 0);
+ monomials(index, 1) = uv(i, 1);
+ monomials(index, 2) = order;
+ }
// other monomials
+ uv = gmshGenerateMonomialsQuadrangle(order + 2);
for (int k = 1; k < order; ++k) {
- fullMatrix<double> uv = gmshGenerateMonomialsQuadrangle(order + 2 - k);
for (int i = 0; i < uv.size1(); ++i, ++index) {
monomials(index, 0) = uv(i, 0);
monomials(index, 1) = uv(i, 1);
@@ -788,10 +779,10 @@ fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
const double p = order + 2;
for (int i = 0; i < points.size1(); ++i) {
- points(i, 2) = points(i, 2) * 1. / p;
- const double duv = -1. + points(i, 2);
- points(i, 0) = duv + points(i, 0) * 2. / p;
- points(i, 1) = duv + points(i, 1) * 2. / p;
+ points(i, 2) = points(i, 2) / p;
+ const double scale = 1. - points(i, 2);
+ points(i, 0) = (-1. + 2. * points(i, 0) / p) * scale;
+ points(i, 1) = (-1. + 2. * points(i, 1) / p) * scale;
}
return points;
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 4b9aeb241729c48280988810ece9f4cf2cec5e6b..c94a3923ba6e3757b9560ef600472562944a6a5e 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -169,7 +169,7 @@ double MetricBasis::getMinR(MElement *el, MetricData *&md, int deg) const
else*/
_computeRmin(*coeff, *jac, RminLag, RminBez, 0, false);
- double betaOpt = beta, minaOpt, maxaOpt = 0., RminBezOpt;
+ double betaOpt = beta, minaOpt = mina, maxaOpt = maxa3, RminBezOpt;
{
/*const */double phi = std::acos(.5*(minK-maxa3*maxa3*maxa3+3*maxa3))/3;
RminBezOpt = (maxa3+2*std::cos(phi+2*M_PI/3))/(maxa3+2*std::cos(phi));
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 9903be53515a43ea52a1cc4a661537ce164e30e8..0fd6fdb956b6d32fdee49fa54b2f6e5b3d29ad33 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -15,182 +15,222 @@
namespace {
// Sub Control Points
- std::vector< fullMatrix<double> > generateSubPointsLine(int order)
- {
- std::vector< fullMatrix<double> > subPoints(2);
+std::vector< fullMatrix<double> > generateSubPointsLine(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(2);
- subPoints[0] = gmshGenerateMonomialsLine(order);
- subPoints[0].scale(.5/order);
+ subPoints[0] = gmshGenerateMonomialsLine(order);
+ subPoints[0].scale(.5/order);
- subPoints[1].copy(subPoints[0]);
- subPoints[1].add(.5);
+ subPoints[1].copy(subPoints[0]);
+ subPoints[1].add(.5);
- return subPoints;
- }
+ return subPoints;
+}
- std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
- {
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
+std::vector< fullMatrix<double> > generateSubPointsTriangle(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
- subPoints[0] = gmshGenerateMonomialsTriangle(order);
- subPoints[0].scale(.5/order);
+ subPoints[0] = gmshGenerateMonomialsTriangle(order);
+ subPoints[0].scale(.5/order);
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
- subPoints[3].copy(subPoints[0]);
- subPoints[3].scale(-1.);
- subPoints[3].add(.5);
+ subPoints[3].copy(subPoints[0]);
+ subPoints[3].scale(-1.);
+ subPoints[3].add(.5);
- return subPoints;
- }
+ return subPoints;
+}
- std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
- {
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
+std::vector< fullMatrix<double> > generateSubPointsQuad(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(4);
+ fullMatrix<double> prox;
- subPoints[0] = gmshGenerateMonomialsQuadrangle(order);
- subPoints[0].scale(.5/order);
+ subPoints[0] = gmshGenerateMonomialsQuadrangle(order);
+ subPoints[0].scale(.5/order);
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
- subPoints[3].copy(subPoints[1]);
- prox.setAsProxy(subPoints[3], 1, 1);
- prox.add(.5);
+ subPoints[3].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[3], 1, 1);
+ prox.add(.5);
- return subPoints;
- }
+ return subPoints;
+}
- std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
- {
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox1;
- fullMatrix<double> prox2;
+std::vector< fullMatrix<double> > generateSubPointsTetrahedron(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox1;
+ fullMatrix<double> prox2;
+
+ subPoints[0] = gmshGenerateMonomialsTetrahedron(order);
+ subPoints[0].scale(.5/order);
+
+ subPoints[1].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[1], 0, 1);
+ prox1.add(.5);
+
+ subPoints[2].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[2], 1, 1);
+ prox1.add(.5);
+
+ subPoints[3].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[3], 2, 1);
+ prox1.add(.5);
+
+ // u := .5-u-w
+ // v := .5-v-w
+ // w := w
+ subPoints[4].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[4], 0, 2);
+ prox1.scale(-1.);
+ prox1.add(.5);
+ prox1.setAsProxy(subPoints[4], 0, 1);
+ prox2.setAsProxy(subPoints[4], 2, 1);
+ prox1.add(prox2, -1.);
+ prox1.setAsProxy(subPoints[4], 1, 1);
+ prox1.add(prox2, -1.);
+
+ // u := u
+ // v := .5-v
+ // w := w+v
+ subPoints[5].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[5], 2, 1);
+ prox2.setAsProxy(subPoints[5], 1, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+ // u := .5-u
+ // v := v
+ // w := w+u
+ subPoints[6].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[6], 2, 1);
+ prox2.setAsProxy(subPoints[6], 0, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+ // u := u+w
+ // v := v+w
+ // w := .5-w
+ subPoints[7].copy(subPoints[0]);
+ prox1.setAsProxy(subPoints[7], 0, 1);
+ prox2.setAsProxy(subPoints[7], 2, 1);
+ prox1.add(prox2);
+ prox1.setAsProxy(subPoints[7], 1, 1);
+ prox1.add(prox2);
+ prox2.scale(-1.);
+ prox2.add(.5);
+
+
+ return subPoints;
+}
- subPoints[0] = gmshGenerateMonomialsTetrahedron(order);
- subPoints[0].scale(.5/order);
+std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox;
- subPoints[1].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[1], 0, 1);
- prox1.add(.5);
+ subPoints[0] = gmshGenerateMonomialsPrism(order);
+ subPoints[0].scale(.5/order);
- subPoints[2].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[2], 1, 1);
- prox1.add(.5);
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
- subPoints[3].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[3], 2, 1);
- prox1.add(.5);
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
- // u := .5-u-w
- // v := .5-v-w
- // w := w
- subPoints[4].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[4], 0, 2);
- prox1.scale(-1.);
- prox1.add(.5);
- prox1.setAsProxy(subPoints[4], 0, 1);
- prox2.setAsProxy(subPoints[4], 2, 1);
- prox1.add(prox2, -1.);
- prox1.setAsProxy(subPoints[4], 1, 1);
- prox1.add(prox2, -1.);
-
- // u := u
- // v := .5-v
- // w := w+v
- subPoints[5].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[5], 2, 1);
- prox2.setAsProxy(subPoints[5], 1, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
- // u := .5-u
- // v := v
- // w := w+u
- subPoints[6].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[6], 2, 1);
- prox2.setAsProxy(subPoints[6], 0, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
-
- // u := u+w
- // v := v+w
- // w := .5-w
- subPoints[7].copy(subPoints[0]);
- prox1.setAsProxy(subPoints[7], 0, 1);
- prox2.setAsProxy(subPoints[7], 2, 1);
- prox1.add(prox2);
- prox1.setAsProxy(subPoints[7], 1, 1);
- prox1.add(prox2);
- prox2.scale(-1.);
- prox2.add(.5);
+ subPoints[3].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[3], 0, 2);
+ prox.scale(-1.);
+ prox.add(.5);
+ subPoints[4].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[4], 2, 1);
+ prox.add(.5);
- return subPoints;
- }
+ subPoints[5].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[5], 2, 1);
+ prox.add(.5);
- std::vector< fullMatrix<double> > generateSubPointsPrism(int order)
- {
- std::vector< fullMatrix<double> > subPoints(8);
- fullMatrix<double> prox;
+ subPoints[6].copy(subPoints[2]);
+ prox.setAsProxy(subPoints[6], 2, 1);
+ prox.add(.5);
- subPoints[0] = gmshGenerateMonomialsPrism(order);
- subPoints[0].scale(.5/order);
+ subPoints[7].copy(subPoints[3]);
+ prox.setAsProxy(subPoints[7], 2, 1);
+ prox.add(.5);
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
+ return subPoints;
+}
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
+std::vector< fullMatrix<double> > generateSubPointsHex(int order)
+{
+ std::vector< fullMatrix<double> > subPoints(8);
+ fullMatrix<double> prox;
- subPoints[3].copy(subPoints[0]);
- prox.setAsProxy(subPoints[3], 0, 2);
- prox.scale(-1.);
- prox.add(.5);
+ subPoints[0] = gmshGenerateMonomialsHexahedron(order);
+ subPoints[0].scale(.5/order);
- subPoints[4].copy(subPoints[0]);
- prox.setAsProxy(subPoints[4], 2, 1);
- prox.add(.5);
+ subPoints[1].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[1], 0, 1);
+ prox.add(.5);
- subPoints[5].copy(subPoints[1]);
- prox.setAsProxy(subPoints[5], 2, 1);
- prox.add(.5);
+ subPoints[2].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[2], 1, 1);
+ prox.add(.5);
- subPoints[6].copy(subPoints[2]);
- prox.setAsProxy(subPoints[6], 2, 1);
- prox.add(.5);
+ subPoints[3].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[3], 1, 1);
+ prox.add(.5);
- subPoints[7].copy(subPoints[3]);
- prox.setAsProxy(subPoints[7], 2, 1);
- prox.add(.5);
+ subPoints[4].copy(subPoints[0]);
+ prox.setAsProxy(subPoints[4], 2, 1);
+ prox.add(.5);
- return subPoints;
- }
+ subPoints[5].copy(subPoints[1]);
+ prox.setAsProxy(subPoints[5], 2, 1);
+ prox.add(.5);
- std::vector< fullMatrix<double> > generateSubPointsHex(int order)
- {
- std::vector< fullMatrix<double> > subPoints(8);
+ subPoints[6].copy(subPoints[2]);
+ prox.setAsProxy(subPoints[6], 2, 1);
+ prox.add(.5);
+
+ subPoints[7].copy(subPoints[3]);
+ prox.setAsProxy(subPoints[7], 2, 1);
+ prox.add(.5);
+
+ return subPoints;
+}
+
+std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
+{
+ if (order == 0) {
+ std::vector< fullMatrix<double> > subPoints(4);
fullMatrix<double> prox;
- subPoints[0] = gmshGenerateMonomialsHexahedron(order);
- subPoints[0].scale(.5/order);
+ subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(0);
+ subPoints[0].scale(.5/(order+2));
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -204,64 +244,14 @@ namespace {
prox.setAsProxy(subPoints[3], 1, 1);
prox.add(.5);
- subPoints[4].copy(subPoints[0]);
- prox.setAsProxy(subPoints[4], 2, 1);
- prox.add(.5);
-
- subPoints[5].copy(subPoints[1]);
- prox.setAsProxy(subPoints[5], 2, 1);
- prox.add(.5);
-
- subPoints[6].copy(subPoints[2]);
- prox.setAsProxy(subPoints[6], 2, 1);
- prox.add(.5);
-
- subPoints[7].copy(subPoints[3]);
- prox.setAsProxy(subPoints[7], 2, 1);
- prox.add(.5);
-
return subPoints;
}
-
- std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
- {
- if (order == 0) {
- std::vector< fullMatrix<double> > subPoints(4);
- fullMatrix<double> prox;
-
- subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
- subPoints[0].scale(.5/(order+2));
-
- subPoints[1].copy(subPoints[0]);
- prox.setAsProxy(subPoints[1], 0, 1);
- prox.add(.5);
-
- subPoints[2].copy(subPoints[0]);
- prox.setAsProxy(subPoints[2], 1, 1);
- prox.add(.5);
-
- subPoints[3].copy(subPoints[1]);
- prox.setAsProxy(subPoints[3], 1, 1);
- prox.add(.5);
-
- return subPoints;
- }
-
- Msg::Error("Subpoints for pyramid p>1 not implemented !");
- return std::vector< fullMatrix<double> >();
-
- /*
+ else {
std::vector< fullMatrix<double> > subPoints(8);
fullMatrix<double> ref, prox;
subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
- int nPts = subPoints[0].size1();
- for (int i = 0; i < nPts; ++i) {
- const double factor = .5 / (order + 2 - subPoints[0](i, 2));
- subPoints[0](i, 0) = subPoints[0](i, 0) * factor;
- subPoints[0](i, 1) = subPoints[0](i, 1) * factor;
- subPoints[0](i, 2) = subPoints[0](i, 2) * .5 / (order + 2);
- }
+ subPoints[0].scale(.5/(order+2));
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -291,6 +281,7 @@ namespace {
prox.setAsProxy(subPoints[7], 2, 1);
prox.add(.5);
+ const int nPts = subPoints[0].size1();
for (int i = 0; i < 8; ++i) {
for (int j = 0; j < nPts; ++j) {
const double factor = (1. - subPoints[i](j, 2));
@@ -300,157 +291,154 @@ namespace {
}
return subPoints;
- */
}
+}
- // Matrices generation
- int nChoosek(int n, int k)
- {
- if (n < k || k < 0) {
- Msg::Error("Wrong argument for combination.");
- return 1;
- }
+// Matrices generation
+int nChoosek(int n, int k)
+{
+ if (n < k || k < 0) {
+ Msg::Error("Wrong argument for combination. (%d, %d)", n, k);
+ return 1;
+ }
- if (k > n/2) k = n-k;
- if (k == 1)
- return n;
- if (k == 0)
- return 1;
+ if (k > n/2) k = n-k;
+ if (k == 1)
+ return n;
+ if (k == 0)
+ return 1;
- int c = 1;
- for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
- return c;
- }
+ int c = 1;
+ for (int i = 1; i <= k; i++, n--) (c *= n) /= i;
+ return c;
+}
- fullMatrix<double> generateBez2LagMatrix
- (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
- int order, int dimSimplex)
- {
- if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
- Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
- exponent.size1(),point.size1(),
- exponent.size2(),point.size2());
- return fullMatrix<double>(1, 1);
- }
+fullMatrix<double> generateBez2LagMatrix
+ (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
+ int order, int dimSimplex)
+{
+ if(exponent.size1() != point.size1() || exponent.size2() != point.size2()){
+ Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
+ exponent.size1(),point.size1(),
+ exponent.size2(),point.size2());
+ return fullMatrix<double>(1, 1);
+ }
- int ndofs = exponent.size1();
- int dim = exponent.size2();
-
- fullMatrix<double> bez2Lag(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- double dd = 1.;
-
- {
- double pointCompl = 1.;
- int exponentCompl = order;
- for (int k = 0; k < dimSimplex; k++) {
- dd *= nChoosek(exponentCompl, (int) exponent(i, k))
- * pow(point(j, k), exponent(i, k));
- pointCompl -= point(j, k);
- exponentCompl -= (int) exponent(i, k);
- }
- dd *= pow(pointCompl, exponentCompl);
+ int ndofs = exponent.size1();
+ int dim = exponent.size2();
+
+ fullMatrix<double> bez2Lag(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double dd = 1.;
+
+ {
+ double pointCompl = 1.;
+ int exponentCompl = order;
+ for (int k = 0; k < dimSimplex; k++) {
+ dd *= nChoosek(exponentCompl, (int) exponent(i, k))
+ * pow(point(j, k), exponent(i, k));
+ pointCompl -= point(j, k);
+ exponentCompl -= (int) exponent(i, k);
}
+ dd *= pow(pointCompl, exponentCompl);
+ }
- for (int k = dimSimplex; k < dim; k++)
- dd *= nChoosek(order, (int) exponent(i, k))
- * pow(point(j, k), exponent(i, k))
- * pow(1. - point(j, k), order - exponent(i, k));
+ for (int k = dimSimplex; k < dim; k++)
+ dd *= nChoosek(order, (int) exponent(i, k))
+ * pow(point(j, k), exponent(i, k))
+ * pow(1. - point(j, k), order - exponent(i, k));
- bez2Lag(j, i) = dd;
- }
+ bez2Lag(j, i) = dd;
}
- return bez2Lag;
}
+ return bez2Lag;
+}
- fullMatrix<double> generateBez2LagMatrixPyramid
- (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
- int order)
- {
- if(exponent.size1() != point.size1() || exponent.size2() != point.size2() ||
- exponent.size2() != 3){
- Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
- exponent.size1(),point.size1(),
- exponent.size2(),point.size2());
- return fullMatrix<double>(1, 1);
- }
+fullMatrix<double> generateBez2LagMatrixPyramid
+ (const fullMatrix<double> &exponent, const fullMatrix<double> &point,
+ int order)
+{
+ if(exponent.size1() != point.size1() || exponent.size2() != point.size2() ||
+ exponent.size2() != 3){
+ Msg::Fatal("Wrong sizes for bez2lag matrix generation %d %d -- %d %d",
+ exponent.size1(),point.size1(),
+ exponent.size2(),point.size2());
+ return fullMatrix<double>(1, 1);
+ }
- int ndofs = exponent.size1();
-
- fullMatrix<double> bez2Lag(ndofs, ndofs);
- for (int i = 0; i < ndofs; i++) {
- for (int j = 0; j < ndofs; j++) {
- bez2Lag(i, j) =
- nChoosek(order, exponent(j, 2))
- * pow_int(point(i, 2), exponent(j, 2))
- * pow_int(1. - point(i, 2), order - exponent(j, 2));
-
- double p = order + 2 - exponent(j, 2);
- double denom = 1. - point(i, 2);
- bez2Lag(i, j) *=
- nChoosek(p, exponent(j, 0))
- * nChoosek(p, exponent(j, 1))
- * pow_int(point(i, 0) / denom, exponent(j, 0))
- * pow_int(point(i, 1) / denom, exponent(j, 1))
- * pow_int(1. - point(i, 0) / denom, p - exponent(j, 0))
- * pow_int(1. - point(i, 1) / denom, p - exponent(j, 1));
- }
+ const int ndofs = exponent.size1();
+
+ fullMatrix<double> bez2Lag(ndofs, ndofs);
+ for (int i = 0; i < ndofs; i++) {
+ for (int j = 0; j < ndofs; j++) {
+ double denom = 1. - point(i, 2);
+ bez2Lag(i, j) =
+ nChoosek(order + 2, exponent(j, 0))
+ * nChoosek(order + 2, exponent(j, 1))
+ * nChoosek(order, exponent(j, 2))
+ * pow_int(point(i, 0) / denom, exponent(j, 0))
+ * pow_int(point(i, 1) / denom, exponent(j, 1))
+ * pow_int(point(i, 2) , exponent(j, 2))
+ * pow_int(1. - point(i, 0) / denom, order + 2 - exponent(j, 0))
+ * pow_int(1. - point(i, 1) / denom, order + 2 - exponent(j, 1))
+ * pow_int(1. - point(i, 2) , order - exponent(j, 2));
}
- return bez2Lag;
}
+ return bez2Lag;
+}
- fullMatrix<double> generateSubDivisor
- (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
- const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
- {
- if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
- Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
- exponents.size1(), lag2Bez.size1(),
- exponents.size1(), lag2Bez.size2());
- return fullMatrix<double>(1, 1);
- }
+fullMatrix<double> generateSubDivisor
+ (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
+ const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
+{
+ if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
+ Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+ exponents.size1(), lag2Bez.size1(),
+ exponents.size1(), lag2Bez.size2());
+ return fullMatrix<double>(1, 1);
+ }
- int nbPts = lag2Bez.size1();
- int nbSubPts = nbPts * subPoints.size();
+ int nbPts = lag2Bez.size1();
+ int nbSubPts = nbPts * subPoints.size();
- fullMatrix<double> intermediate2(nbPts, nbPts);
- fullMatrix<double> subDivisor(nbSubPts, nbPts);
+ fullMatrix<double> intermediate2(nbPts, nbPts);
+ fullMatrix<double> subDivisor(nbSubPts, nbPts);
- for (unsigned int i = 0; i < subPoints.size(); i++) {
- fullMatrix<double> intermediate1 =
- generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
- lag2Bez.mult(intermediate1, intermediate2);
- subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
- }
- return subDivisor;
+ for (unsigned int i = 0; i < subPoints.size(); i++) {
+ fullMatrix<double> intermediate1 =
+ generateBez2LagMatrix(exponents, subPoints[i], order, dimSimplex);
+ lag2Bez.mult(intermediate1, intermediate2);
+ subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
}
+ return subDivisor;
+}
- fullMatrix<double> generateSubDivisorPyramid
- (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
- const fullMatrix<double> &lag2Bez, int order)
- {
- if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
- Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
- exponents.size1(), lag2Bez.size1(),
- exponents.size1(), lag2Bez.size2());
- return fullMatrix<double>(1, 1);
- }
+fullMatrix<double> generateSubDivisorPyramid
+ (const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
+ const fullMatrix<double> &lag2Bez, int order)
+{
+ if (exponents.size1() != lag2Bez.size1() || exponents.size1() != lag2Bez.size2()){
+ Msg::Fatal("Wrong sizes for Bezier Divisor %d %d -- %d %d",
+ exponents.size1(), lag2Bez.size1(),
+ exponents.size1(), lag2Bez.size2());
+ return fullMatrix<double>(1, 1);
+ }
- int nbPts = lag2Bez.size1();
- int nbSubPts = nbPts * subPoints.size();
+ int nbPts = lag2Bez.size1();
+ int nbSubPts = nbPts * subPoints.size();
- fullMatrix<double> intermediate2(nbPts, nbPts);
- fullMatrix<double> subDivisor(nbSubPts, nbPts);
+ fullMatrix<double> intermediate2(nbPts, nbPts);
+ fullMatrix<double> subDivisor(nbSubPts, nbPts);
- for (unsigned int i = 0; i < subPoints.size(); i++) {
- fullMatrix<double> intermediate1 =
- generateBez2LagMatrixPyramid(exponents, subPoints[i], order);
- lag2Bez.mult(intermediate1, intermediate2);
- subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
- }
- return subDivisor;
+ for (unsigned int i = 0; i < subPoints.size(); i++) {
+ fullMatrix<double> intermediate1 =
+ generateBez2LagMatrixPyramid(exponents, subPoints[i], order);
+ lag2Bez.mult(intermediate1, intermediate2);
+ subDivisor.copy(intermediate2, 0, nbPts, 0, nbPts, i*nbPts, 0);
}
+ return subDivisor;
+}
}
void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
@@ -531,23 +519,23 @@ void bezierBasis::_construct(int parentType, int p)
dim = 3;
numLagCoeff = 8;
_exponents = JacobianBasis::generateJacMonomialsPyramid(order);
- if (order < 2) {
- subPoints = generateSubPointsPyr(order);
- numDivisions = static_cast<int>(subPoints.size());
- }
- else {
- numDivisions = 8;
- static int a = 0;
- if (++a == 1) Msg::Warning("Subdivision not available for pyramid p>1");
- }
- fullMatrix<double> bezierPoints = _exponents;
- bezierPoints.scale(1. / (order + 2));
+ subPoints = generateSubPointsPyr(order);
+ numDivisions = static_cast<int>(subPoints.size());
+
+ fullMatrix<double> bezierPoints;
+ bezierPoints.resize(_exponents.size1(), _exponents.size2());
+ const double p = order + 2;
+ for (int i = 0; i < bezierPoints.size1(); ++i) {
+ bezierPoints(i, 2) = _exponents(i, 2) / p;
+ const double scale = 1. - bezierPoints(i, 2);
+ bezierPoints(i, 0) = _exponents(i, 0) / p * scale;
+ bezierPoints(i, 1) = _exponents(i, 1) / p * scale;
+ }
matrixBez2Lag = generateBez2LagMatrixPyramid(_exponents, bezierPoints, order);
matrixBez2Lag.invert(matrixLag2Bez);
- if (order < 2)
- subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez, order);
+ subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez, order);
return;
}
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index c71682284a3412358353509df2e08af4515e039f..657ccbc53a1491efdd729a22f88e4f0a96e70823 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -110,7 +110,9 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
"\n"
"- Recompute = {0,1}: If the mesh has changed, set to 1 to recompute the bounds.\n"
"\n"
- "- Tolerance = ]0, 1[: Tolerance on the computation of min(R) or min(J)";
+ "- Tolerance = ]0, 1[: Tolerance on the computation of min(R) or min(J). "
+ "It should be at most 0.01 but it can be set to 1 to just check the validity of "
+ "the mesh.";
}
// Execution