diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 95656b01a429afd4813eb07de7aedd295cc74673..6f1584e1bdaedafe494bef6c8994d53c6b296978 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -19,7 +19,8 @@ JacobianBasis::JacobianBasis(int tag)
const int parentType = ElementType::ParentTypeFromTag(tag);
const int order = ElementType::OrderFromTag(tag);
- int jacType = 0, primJacType = 0;
+ //int jacType = 0,
+ int primJacType = 0;
/*if (parentType == TYPE_PYR) {
switch (tag) {
@@ -33,31 +34,59 @@ JacobianBasis::JacobianBasis(int tag)
}
else {*/
int jacobianOrder, primJacobianOrder;
+
switch (parentType) {
- case TYPE_PNT : jacobianOrder = 0; primJacobianOrder = 0; break;
- case TYPE_LIN : jacobianOrder = order - 1; primJacobianOrder = 0; break;
- case TYPE_TRI : jacobianOrder = 2 * (order - 1); primJacobianOrder = 0; break;
- case TYPE_QUA : jacobianOrder = 2 * order - 1; primJacobianOrder = 1; break;
- case TYPE_TET : jacobianOrder = 3 * (order - 1); primJacobianOrder = 0; break;
- case TYPE_PRI : jacobianOrder = 3 * order - 1; primJacobianOrder = 2; break;
- case TYPE_HEX : jacobianOrder = 3 * order - 1; primJacobianOrder = 2; break;
+ case TYPE_PNT :
+ primJacobianOrder = 0;
+ jacobianOrder = 0;
+ lagPoints = gmshGeneratePointsLine(0);
+ break;
+ case TYPE_LIN :
+ primJacobianOrder = 0;
+ jacobianOrder = order - 1;
+ lagPoints = gmshGeneratePointsLine(order);
+ break;
+ case TYPE_TRI :
+ primJacobianOrder = 0;
+ jacobianOrder = 2*order - 2;
+ lagPoints = gmshGeneratePointsTriangle(order);
+ break;
+ case TYPE_QUA :
+ primJacobianOrder = 1;
+ jacobianOrder = 2*order - 1;
+ lagPoints = gmshGeneratePointsQuadrangle(order);
+ break;
+ case TYPE_TET :
+ primJacobianOrder = 0;
+ jacobianOrder = 3*order - 3;
+ lagPoints = gmshGeneratePointsTetrahedron(order);
+ break;
+ case TYPE_PRI :
+ primJacobianOrder = 2;
+ jacobianOrder = 3*order - 1;
+ lagPoints = gmshGeneratePointsPrism(order);
+ break;
+ case TYPE_HEX :
+ primJacobianOrder = 2;
+ jacobianOrder = 3*order - 1;
+ lagPoints = gmshGeneratePointsHexahedron(order);
+ break;
default :
Msg::Error("Unknown Jacobian function space for element type %d", tag);
jacobianOrder = 0;
break;
}
- jacType = ElementType::getTag(parentType, jacobianOrder, false);
primJacType = ElementType::getTag(parentType, primJacobianOrder, false);
//}
// Store Bezier basis
- bezier = BasisFactory::getBezierBasis(jacType);
+ bezier = BasisFactory::getBezierBasis(parentType, jacobianOrder);
// Store shape function gradients of mapping at Jacobian nodes
const nodalBasis *mapBasis = BasisFactory::getNodalBasis(tag);
fullMatrix<double> allDPsi;
- mapBasis->df(getPoints(), allDPsi);
- numJacNodes = getPoints().size1();
+ mapBasis->df(lagPoints, allDPsi);
+ numJacNodes = lagPoints.size1();
numMapNodes = allDPsi.size1();
gradShapeMatX.resize(numJacNodes,numMapNodes);
gradShapeMatY.resize(numJacNodes,numMapNodes);
@@ -73,7 +102,7 @@ JacobianBasis::JacobianBasis(int tag)
const nodalBasis *primJacBasis = BasisFactory::getNodalBasis(primJacType);
numPrimJacNodes = primJacBasis->getNumShapeFunctions();
matrixPrimJac2Jac.resize(numJacNodes,numPrimJacNodes);
- primJacBasis->f(getPoints(),matrixPrimJac2Jac);
+ primJacBasis->f(lagPoints, matrixPrimJac2Jac); //FIXME not too much sampling points ??
// Compute shape function gradients of primary mapping at barycenter,
// in order to compute normal to straight element
@@ -193,7 +222,7 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
- switch (bezier->dim) {
+ switch (bezier->getDim()) {
case 1 : {
fullMatrix<double> normals(2,3);
@@ -225,7 +254,7 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVe
void JacobianBasis::getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
- switch (bezier->dim) {
+ switch (bezier->getDim()) {
case 1 : {
fullMatrix<double> normals(2,3);
@@ -262,7 +291,7 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &normals, fullVector<double> &jacobian) const
{
- switch (bezier->dim) {
+ switch (bezier->getDim()) {
case 0 : {
for (int i = 0; i < numJacNodes; i++) jacobian(i) = 1.;
@@ -270,7 +299,7 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ,
}
case 1 : {
- fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3);
+ fullMatrix<double> dxyzdX(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
for (int i = 0; i < numJacNodes; i++) {
const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
@@ -319,7 +348,7 @@ void JacobianBasis::getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
fullMatrix<double> &JDJ) const
{
- switch (bezier->dim) {
+ switch (bezier->getDim()) {
case 0 : {
for (int i = 0; i < numJacNodes; i++) {
@@ -473,7 +502,7 @@ void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesX, const fu
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
- switch (bezier->dim) {
+ switch (bezier->getDim()) {
case 0 : {
const int numEl = nodesX.size2();
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index d5f34792bfe5d0fe36e62ba3b22cfead7b63729c..b38e78aa73dba7b8e19f1a5442a4c977172c3eb5 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -15,11 +15,11 @@
class JacobianBasis {
private:
const bezierBasis *bezier;
+ fullMatrix<double> lagPoints; // sampling point
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
fullVector<double> primGradShapeBarX, primGradShapeBarY, primGradShapeBarZ;
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numMapNodes, numPrimJacNodes, numPrimMapNodes;
- inline const fullMatrix<double> &getPoints() const { return bezier->lagPoints; }
public :
JacobianBasis(int tag);
@@ -29,9 +29,9 @@ class JacobianBasis {
inline int getNumMapNodes() const { return numMapNodes; }
inline int getNumPrimJacNodes() const { return numPrimJacNodes; }
inline int getNumPrimMapNodes() const { return numPrimMapNodes; }
- inline int getNumDivisions() const { return bezier->numDivisions; }
+ inline int getNumDivisions() const { return bezier->getNumDivision(); }
inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
- inline int getNumLagPts() const { return bezier->numLagPts; }
+ inline int getNumLagCoeff() const { return bezier->getNumLagCoeff(); }
//
double getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index f8884abebae5ff040d7676ba00167dea36d8d742..7db0b0904cc783aa0cc29cc91fb546d17c623730 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -336,62 +336,55 @@ void bezierBasis::_construct(int parentType, int p)
switch (parentType) {
case TYPE_PNT :
dim = 0;
- numLagPts = 1;
+ numLagCoeff = 1;
dimSimplex = 0;
exponents = gmshGenerateMonomialsLine(0);
- lagPoints = gmshGeneratePointsLine(0);
break;
case TYPE_LIN : {
dim = 1;
- numLagPts = 2;
+ numLagCoeff = 2;
dimSimplex = 0;
exponents = gmshGenerateMonomialsLine(order);
- lagPoints = gmshGeneratePointsLine(order);
subPoints = generateSubPointsLine(order);
break;
}
case TYPE_TRI : {
dim = 2;
- numLagPts = 3;
+ numLagCoeff = 3;
dimSimplex = 2;
exponents = gmshGenerateMonomialsTriangle(order);
- lagPoints = gmshGeneratePointsTriangle(order);
subPoints = generateSubPointsTriangle(order);
break;
}
case TYPE_QUA : {
dim = 2;
- numLagPts = 4;
+ numLagCoeff = 4;
dimSimplex = 0;
exponents = gmshGenerateMonomialsQuadrangle(order);
- lagPoints = gmshGeneratePointsQuadrangle(order);
subPoints = generateSubPointsQuad(order);
break;
}
case TYPE_TET : {
dim = 3;
- numLagPts = 4;
+ numLagCoeff = 4;
dimSimplex = 3;
exponents = gmshGenerateMonomialsTetrahedron(order);
- lagPoints = gmshGeneratePointsTetrahedron(order);
subPoints = generateSubPointsTetrahedron(order);
break;
}
case TYPE_PRI : {
dim = 3;
- numLagPts = 6;
+ numLagCoeff = 6;
dimSimplex = 2;
exponents = gmshGenerateMonomialsPrism(order);
- lagPoints = gmshGeneratePointsPrism(order);
subPoints = generateSubPointsPrism(order);
break;
}
case TYPE_HEX : {
dim = 3;
- numLagPts = 8;
+ numLagCoeff = 8;
dimSimplex = 0;
exponents = gmshGenerateMonomialsHexahedron(order);
- lagPoints = gmshGeneratePointsHexahedron(order);
subPoints = generateSubPointsHex(order);
break;
}
@@ -400,10 +393,9 @@ void bezierBasis::_construct(int parentType, int p)
"reverting to TET_1", parentType);
dim = 3;
order = 0;
- numLagPts = 4;
+ numLagCoeff = 4;
dimSimplex = 3;
exponents = gmshGenerateMonomialsTetrahedron(order);
- lagPoints = gmshGeneratePointsTetrahedron(order);
subPoints = generateSubPointsTetrahedron(order);
break;
}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index 62aca79d0810d1c243f81c3ec0f832a119d1bab8..4e0c519383de1637d4aa73da22d8ad7d06e0e44b 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -14,18 +14,18 @@
class MElement;
class bezierBasis {
- public:
- //private :
- // the 'numLagPts' first exponents are related to 'Lagrangian' values
+ private :
+ // the 'numLagCoeff' first exponents are related to 'Lagrangian' values
+ int numLagCoeff;
int dim, order;
- int numLagPts;
int numDivisions;
- fullMatrix<double> subDivisor;
public :
- fullMatrix<double> lagPoints; // sampling point
- fullMatrix<double> matrixLag2Bez, matrixBez2Lag;
+ fullMatrix<double> matrixLag2Bez;
+ fullMatrix<double> matrixBez2Lag;
+ fullMatrix<double> subDivisor;
+ // Constructors
inline bezierBasis(int tag) {
_construct(ElementType::ParentTypeFromTag(tag), ElementType::OrderFromTag(tag));
}
@@ -34,10 +34,10 @@ class bezierBasis {
}
// get methods
- inline int getDim() {return dim;}
- inline int getOrder() {return order;}
- inline int getNumLagPts() {return numLagPts;}
- inline int getNumDivision() {return numDivisions;}
+ inline int getDim() const {return dim;}
+ inline int getOrder() const {return order;}
+ inline int getNumLagCoeff() const {return numLagCoeff;}
+ inline int getNumDivision() const {return numDivisions;}
private :
void _construct(int parendtType, int order);
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 1e4f5438ce8d92dc3a847a0419e1060c50fc2704..182e76d5dabcd734618d2e971d7be0b86efcef48 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -339,7 +339,7 @@ int GMSH_AnalyseCurvedMeshPlugin::subDivision(const JacobianBasis *jb,
jb->subDivisor(jacobian, newJacobian);
for (int i = 0; i < jb->getNumDivisions(); i++)
- for (int j = 0; j < jb->getNumLagPts(); j++)
+ for (int j = 0; j < jb->getNumLagCoeff(); j++)
if (newJacobian(i * jb->getNumJacNodes() + j) <= _jacBreak)
return -1;
@@ -727,7 +727,7 @@ BezierJacobian::BezierJacobian(fullVector<double> &v, const JacobianBasis *jfs,
_minJ = _maxJ = v(0);
int i = 1;
- for (; i < jfs->getNumLagPts(); i++) {
+ for (; i < jfs->getNumLagCoeff(); i++) {
if (_minJ > v(i)) _minJ = v(i);
if (_maxJ < v(i)) _maxJ = v(i);
}