diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index a58d3fbc472c971a4d627b4a278fb415def2363d..8975c7006dfcfda5547138548ba55d7db588cebb 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -18,53 +18,37 @@ JacobianBasis::JacobianBasis(int tag)
{
const int parentType = ElementType::ParentTypeFromTag(tag);
const int order = ElementType::OrderFromTag(tag);
+ const int jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
+ const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
- int jacobianOrder, primJacobianOrder;
switch (parentType) {
case TYPE_PNT :
- primJacobianOrder = 0;
- jacobianOrder = 0;
lagPoints = gmshGeneratePointsLine(0);
break;
case TYPE_LIN :
- primJacobianOrder = 0;
- jacobianOrder = order - 1;
lagPoints = gmshGeneratePointsLine(jacobianOrder);
break;
case TYPE_TRI :
- primJacobianOrder = 0;
- jacobianOrder = 2*order - 2;
lagPoints = gmshGeneratePointsTriangle(jacobianOrder);
break;
case TYPE_QUA :
- primJacobianOrder = 1;
- jacobianOrder = 2*order - 1;
lagPoints = gmshGeneratePointsQuadrangle(jacobianOrder);
break;
case TYPE_TET :
- primJacobianOrder = 0;
- jacobianOrder = 3*order - 3;
lagPoints = gmshGeneratePointsTetrahedron(jacobianOrder);
break;
case TYPE_PRI :
- primJacobianOrder = 2;
- jacobianOrder = 3*order - 1;
lagPoints = gmshGeneratePointsPrism(jacobianOrder);
break;
case TYPE_HEX :
- primJacobianOrder = 2;
- jacobianOrder = 3*order - 1;
lagPoints = gmshGeneratePointsHexahedron(jacobianOrder);
break;
case TYPE_PYR :
- primJacobianOrder = 0;
- jacobianOrder = 3*order - 3;
lagPoints = generateJacPointsPyramid(jacobianOrder);
break;
default :
Msg::Error("Unknown Jacobian function space for element type %d", tag);
- jacobianOrder = 0;
- break;
+ return;
}
numJacNodes = lagPoints.size1();
@@ -680,3 +664,36 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
}
return monomials;
}
+
+fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
+{
+ fullMatrix<double> points = generateJacMonomialsPyramid(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;
+ }
+
+ return points;
+}
+
+int JacobianBasis::jacobianOrder(int parentType, int order) {
+ switch (parentType) {
+ case TYPE_PNT : return 0;
+ case TYPE_LIN : return order - 1;
+ case TYPE_TRI : return 2*order - 2;
+ case TYPE_QUA : return 2*order - 1;
+ case TYPE_TET : return 3*order - 3;
+ case TYPE_PRI : return 3*order - 1;
+ case TYPE_HEX : return 3*order - 1;
+ case TYPE_PYR : return 3*order - 3;
+ // note : for the pyramid, the jacobian space is
+ // different from the space of the mapping
+ default :
+ Msg::Error("Unknown element type %d, return order 0", parentType);
+ return 0;
+ }
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index df8ebac767eb7df3b1e782d8d7fa8cd6f8664c30..5aa77258dca7432f05553b2177dd22e9af5745ed 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -67,15 +67,9 @@ class JacobianBasis {
}
//
+ static int jacobianOrder(int parentType, int order);
static fullMatrix<double> generateJacMonomialsPyramid(int order);
- static inline fullMatrix<double> generateJacPointsPyramid(int order) {
- fullMatrix<double> points = generateJacMonomialsPyramid(order);
-
- if (order == 0) return points;
-
- points.scale(1. / (order+2.));
- return points;
- }
+ static fullMatrix<double> generateJacPointsPyramid(int order);
};
#endif
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 7db0b0904cc783aa0cc29cc91fb546d17c623730..cf2d14c8e413d63e9b30f4c83ba82de74c950fc6 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -13,7 +13,6 @@
#include "BasisFactory.h"
#include "Numeric.h"
-
// Sub Control Points
static std::vector< fullMatrix<double> > generateSubPointsLine(int order)
{
@@ -223,6 +222,81 @@ static std::vector< fullMatrix<double> > generateSubPointsHex(int order)
return subPoints;
}
+static 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;
+ }
+
+ 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[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);
+
+ 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);
+
+ for (int i = 0; i < 8; ++i) {
+ for (int j = 0; j < nPts; ++j) {
+ const double factor = (1. - subPoints[i](j, 2));
+ subPoints[i](j, 0) = subPoints[i](j, 0) * factor;
+ subPoints[i](j, 1) = subPoints[i](j, 1) * factor;
+ }
+ }
+
+ return subPoints;
+}
+
// Matrices generation
static int nChoosek(int n, int k)
{
@@ -246,9 +320,8 @@ static 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 Bezier coefficients generation %d %d -- %d %d",
+ 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);
@@ -285,6 +358,42 @@ static fullMatrix<double> generateBez2LagMatrix
return bez2Lag;
}
+static 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(point(i, 2), exponent(j, 2)) //FIXME use pow_int
+ * pow(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(point(i, 0) / denom, exponent(j, 0))
+ * pow(point(i, 1) / denom, exponent(j, 1))
+ * pow(1. - point(i, 0) / denom, p - exponent(j, 0))
+ * pow(1. - point(i, 1) / denom, p - exponent(j, 1));
+ }
+ }
+ return bez2Lag;
+}
+
static fullMatrix<double> generateSubDivisor
(const fullMatrix<double> &exponents, const std::vector< fullMatrix<double> > &subPoints,
const fullMatrix<double> &lag2Bez, int order, int dimSimplex)
@@ -311,102 +420,130 @@ static fullMatrix<double> generateSubDivisor
return subDivisor;
}
+static 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();
+
+ 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;
+}
+
void bezierBasis::_construct(int parentType, int p)
{
order = p;
+ fullMatrix<double> exponents;
+ std::vector< fullMatrix<double> > subPoints;
- /*if (parentType == TYPE_PYR) {
+ if (parentType == TYPE_PYR) {
dim = 3;
- numLagPts = 5;
- bezierPoints = gmshGeneratePointsPyramid(order,false);
- lagPoints = bezierPoints;
- matrixLag2Bez.resize(bezierPoints.size1(), bezierPoints.size1(), false);
- matrixBez2Lag.resize(bezierPoints.size1(), bezierPoints.size1(), false);
- for (int i=0; i<matrixBez2Lag.size1(); ++i) {
- matrixLag2Bez(i,i) = 1.;
- matrixBez2Lag(i,i) = 1.;
- }
- // TODO: Implement subdidivsor
- }
- else {*/
- int dimSimplex;
- fullMatrix<double> exponents;
- std::vector< fullMatrix<double> > subPoints;
-
- switch (parentType) {
- case TYPE_PNT :
- dim = 0;
- numLagCoeff = 1;
- dimSimplex = 0;
- exponents = gmshGenerateMonomialsLine(0);
- break;
- case TYPE_LIN : {
- dim = 1;
- numLagCoeff = 2;
- dimSimplex = 0;
- exponents = gmshGenerateMonomialsLine(order);
- subPoints = generateSubPointsLine(order);
- break;
- }
- case TYPE_TRI : {
- dim = 2;
- numLagCoeff = 3;
- dimSimplex = 2;
- exponents = gmshGenerateMonomialsTriangle(order);
- subPoints = generateSubPointsTriangle(order);
- break;
- }
- case TYPE_QUA : {
- dim = 2;
- numLagCoeff = 4;
- dimSimplex = 0;
- exponents = gmshGenerateMonomialsQuadrangle(order);
- subPoints = generateSubPointsQuad(order);
- break;
- }
- case TYPE_TET : {
- dim = 3;
- numLagCoeff = 4;
- dimSimplex = 3;
- exponents = gmshGenerateMonomialsTetrahedron(order);
- subPoints = generateSubPointsTetrahedron(order);
- break;
- }
- case TYPE_PRI : {
- dim = 3;
- numLagCoeff = 6;
- dimSimplex = 2;
- exponents = gmshGenerateMonomialsPrism(order);
- subPoints = generateSubPointsPrism(order);
- break;
- }
- case TYPE_HEX : {
- dim = 3;
- numLagCoeff = 8;
- dimSimplex = 0;
- exponents = gmshGenerateMonomialsHexahedron(order);
- subPoints = generateSubPointsHex(order);
- break;
- }
- default : {
- Msg::Error("Unknown function space of parentType %d : "
- "reverting to TET_1", parentType);
- dim = 3;
- order = 0;
- numLagCoeff = 4;
- dimSimplex = 3;
- exponents = gmshGenerateMonomialsTetrahedron(order);
- subPoints = generateSubPointsTetrahedron(order);
- break;
- }
- }
+ numLagCoeff = 8;
+ exponents = JacobianBasis::generateJacMonomialsPyramid(order);
+ subPoints = generateSubPointsPyr(order);
+
+ numDivisions = static_cast<int>(subPoints.size());
fullMatrix<double> bezierPoints = exponents;
- bezierPoints.scale(1./order);
+ bezierPoints.scale(1. / (order + 2));
- numDivisions = static_cast<int>(subPoints.size());
- matrixBez2Lag = generateBez2LagMatrix(exponents, bezierPoints, order, dimSimplex);
+ matrixBez2Lag = generateBez2LagMatrixPyramid(exponents, bezierPoints, order);
matrixBez2Lag.invert(matrixLag2Bez);
- subDivisor = generateSubDivisor(exponents, subPoints, matrixLag2Bez, order, dimSimplex);
- //}
+ subDivisor = generateSubDivisorPyramid(exponents, subPoints, matrixLag2Bez, order);
+ return;
+ }
+
+ int dimSimplex;
+ switch (parentType) {
+ case TYPE_PNT :
+ dim = 0;
+ numLagCoeff = 1;
+ dimSimplex = 0;
+ exponents = gmshGenerateMonomialsLine(0);
+ subPoints.push_back(gmshGeneratePointsLine(0));
+ break;
+ case TYPE_LIN : {
+ dim = 1;
+ numLagCoeff = 2;
+ dimSimplex = 0;
+ exponents = gmshGenerateMonomialsLine(order);
+ subPoints = generateSubPointsLine(order);
+ break;
+ }
+ case TYPE_TRI : {
+ dim = 2;
+ numLagCoeff = 3;
+ dimSimplex = 2;
+ exponents = gmshGenerateMonomialsTriangle(order);
+ subPoints = generateSubPointsTriangle(order);
+ break;
+ }
+ case TYPE_QUA : {
+ dim = 2;
+ numLagCoeff = 4;
+ dimSimplex = 0;
+ exponents = gmshGenerateMonomialsQuadrangle(order);
+ subPoints = generateSubPointsQuad(order);
+ break;
+ }
+ case TYPE_TET : {
+ dim = 3;
+ numLagCoeff = 4;
+ dimSimplex = 3;
+ exponents = gmshGenerateMonomialsTetrahedron(order);
+ subPoints = generateSubPointsTetrahedron(order);
+ break;
+ }
+ case TYPE_PRI : {
+ dim = 3;
+ numLagCoeff = 6;
+ dimSimplex = 2;
+ exponents = gmshGenerateMonomialsPrism(order);
+ subPoints = generateSubPointsPrism(order);
+ break;
+ }
+ case TYPE_HEX : {
+ dim = 3;
+ numLagCoeff = 8;
+ dimSimplex = 0;
+ exponents = gmshGenerateMonomialsHexahedron(order);
+ subPoints = generateSubPointsHex(order);
+ break;
+ }
+ default : {
+ Msg::Error("Unknown function space of parentType %d : "
+ "reverting to TET_1", parentType);
+ dim = 3;
+ order = 0;
+ numLagCoeff = 4;
+ dimSimplex = 3;
+ exponents = gmshGenerateMonomialsTetrahedron(order);
+ subPoints = generateSubPointsTetrahedron(order);
+ break;
+ }
+ }
+ numDivisions = static_cast<int>(subPoints.size());
+
+ fullMatrix<double> bezierPoints = exponents;
+ bezierPoints.scale(1./order);
+
+ matrixBez2Lag = generateBez2LagMatrix(exponents, bezierPoints, order, dimSimplex);
+ matrixBez2Lag.invert(matrixLag2Bez);
+ subDivisor = generateSubDivisor(exponents, subPoints, matrixLag2Bez, order, dimSimplex);
}