diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index b07ae2fdc2f4a4dae308232a098a3b8e7dd84fa8..ab0e1cb43af5f31ccd0b84cd4d5c838dfdff5d81 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -113,21 +113,20 @@ void MElement::scaledJacRange(double &jmin, double &jmax)
jmin = jmax = 1.0;
#if defined(HAVE_MESH)
if (getPolynomialOrder() == 1) return;
- const JacobianBasis *jac = getJacobianFuncSpace(); // Calc signed Jacobian
- const bezierBasis *bez = jac->bezier;
- const int numJacNodes = bez->points.size1();
- fullVector<double> Ji(numJacNodes);
- jac->getSignedJacobian(this,Ji);
- const JacobianBasis *jac1 = getJacobianFuncSpace(1); // Calc signed primary Jacobian
- const int numJac1Nodes = jac1->bezier->points.size1();
- fullVector<double> J1i(numJac1Nodes), J1Ji(numJacNodes);
- jac1->getSignedJacobian(this,J1i);
- for (int i=0; i<numJac1Nodes; i++) if (J1i(i) == 0) std::cout << "DBGTT: el " << getNum() << ", Ji1(" << i << ") == 0!\n";
- jac->primJac2Jac.mult(J1i,J1Ji);
- fullVector<double> SJi(numJacNodes); // Calc scaled Jacobian
+ const JacobianBasis *jac = getJacobianFuncSpace(),*jac1 = getJacobianFuncSpace(1); // Jac. and prim. Jac. basis
+ const int numJacNodes = jac->getNumJacNodes(), numJac1Nodes = jac1->getNumJacNodes();
+ const int numMapNodes = jac->getNumMapNodes(), numMap1Nodes = jac1->getNumMapNodes();
+ const int dim = getDim();
+ fullMatrix<double> nodesXYZ(numMapNodes,dim), nodesXYZ1(numMap1Nodes,dim);
+ getNodesCoord(nodesXYZ);
+ nodesXYZ1.copy(nodesXYZ,0,numMap1Nodes,0,dim,0,0);
+ fullVector<double> Ji(numJacNodes), J1i(numJac1Nodes), J1Ji(numJacNodes);
+ jac->getSignedJacobian(nodesXYZ,Ji);
+ jac1->getSignedJacobian(nodesXYZ1,J1i);
+ jac->primJac2Jac(J1i,J1Ji);
+ fullVector<double> SJi(numJacNodes), Bi(numJacNodes); // Calc scaled Jacobian -> Bezier
for (int i=0; i<numJacNodes; i++) SJi(i) = Ji(i)/J1Ji(i);
- fullVector<double> Bi(numJacNodes); // Transform to Bezier basis
- bez->matrixLag2Bez.mult(SJi,Bi);
+ jac->lag2Bez(SJi,Bi);
jmin = *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
jmax = *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
#endif
@@ -412,6 +411,34 @@ double MElement::getPrimaryJacobian(double u, double v, double w, double jac[3][
return _computeDeterminantAndRegularize(this, jac);
}
+void MElement::getSignedJacobian(fullVector<double> &jacobian, int o)
+{
+ const int dim = getDim(), numNodes = getNumVertices();
+ fullMatrix<double> nodesXYZ(numNodes,dim);
+ getNodesCoord(nodesXYZ);
+ getJacobianFuncSpace(o)->getSignedJacobian(nodesXYZ,jacobian);
+}
+
+void MElement::getNodesCoord(fullMatrix<double> &nodesXYZ)
+{
+ const int dim = getDim(), numNodes = getNumShapeFunctions();
+ if (dim <= 1)
+ for (int i = 0; i < numNodes; i++) nodesXYZ(i,0) = getShapeFunctionNode(i)->x();
+ else if (dim == 2)
+ for (int i = 0; i < numNodes; i++) {
+ MVertex *v = getShapeFunctionNode(i);
+ nodesXYZ(i,0) = v->x();
+ nodesXYZ(i,1) = v->y();
+ }
+ else if (dim == 3)
+ for (int i = 0; i < numNodes; i++) {
+ MVertex *v = getShapeFunctionNode(i);
+ nodesXYZ(i,0) = v->x();
+ nodesXYZ(i,1) = v->y();
+ nodesXYZ(i,2) = v->z();
+ }
+}
+
void MElement::pnt(double u, double v, double w, SPoint3 &p)
{
double x = 0., y = 0., z = 0.;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index c6f3a61af706a66e862f80d6dfb5f43d0ba5673b..1b56081c8cc0272dc219b1f8fd6b3285bdf82ba4 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -264,6 +264,8 @@ class MElement
{
double jac[3][3]; return getJacobian(u, v, w, jac);
}
+ void getSignedJacobian(fullVector<double> &jacobian, int o = -1);
+ void getNodesCoord(fullMatrix<double> &nodesXYZ);
virtual int getNumShapeFunctions(){ return getNumVertices(); }
virtual int getNumPrimaryShapeFunctions(){ return getNumPrimaryVertices(); }
virtual MVertex *getShapeFunctionNode(int i){ return getVertex(i); }
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 6265eb91fe33eab526a17eaaf100ec0350d1803b..23cf8e8e6c3cfcbe52995260c202bc2e609ff9cd 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -8,6 +8,7 @@
#include "Numeric.h"
#include "Context.h"
#include "BasisFactory.h"
+#include "JacobianBasis.h"
#if defined(HAVE_MESH)
#include "qualityMeasures.h"
@@ -51,13 +52,13 @@ void MTetrahedron10::scaledJacRange(double &jmin, double &jmax){
const int nbBezT = 20;
static int done = 0;
static fullMatrix<double> gsf[nbBezT];
- const bezierBasis *jac_basis = getJacobianFuncSpace()->bezier;
+ const JacobianBasis *jac_basis = getJacobianFuncSpace();
if (!done){
double gs[nbNodT][3];
- for (int i=0;i<jac_basis->points.size1();i++){
- const double u = jac_basis->points(i,0);
- const double v = jac_basis->points(i,1);
- const double w = jac_basis->points(i,2);
+ for (int i=0;i<jac_basis->getNumJacNodes();i++){
+ const double u = jac_basis->getPoints()(i,0);
+ const double v = jac_basis->getPoints()(i,1);
+ const double w = jac_basis->getPoints()(i,2);
getGradShapeFunctions(u,v,w,gs);
fullMatrix<double> a (nbNodT,3);
for (int j=0;j < nbNodT;j++){
@@ -104,7 +105,7 @@ void MTetrahedron10::scaledJacRange(double &jmin, double &jmax){
Ji(i) = dJ * ss;
}
static fullVector<double> Bi( nbBezT );
- jac_basis->matrixLag2Bez.mult(Ji,Bi);
+ jac_basis->lag2Bez(Ji,Bi);
// printf("%22.15E\n",*std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size()));
jmin= *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
jmax= *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 37b5f13bf460590ebaa3174e5a8a09a35b1ff3f6..46a8a908c283299c3f9e209a537a0e36da716138 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -337,11 +337,11 @@ double mesh_functional_distorsion_p2_exact(MTriangle *t)
double mesh_functional_distorsion_pN(MElement *t)
{
- const bezierBasis *jac = t->getJacobianFuncSpace()->bezier;
- fullVector<double>Ji(jac->points.size1());
- for (int i=0;i<jac->points.size1();i++){
- double u = jac->points(i,0);
- double v = jac->points(i,1);
+ const JacobianBasis *jac = t->getJacobianFuncSpace();
+ fullVector<double>Ji(jac->getNumJacNodes());
+ for (int i=0;i<jac->getNumJacNodes();i++){
+ double u = jac->getPoints()(i,0);
+ double v = jac->getPoints()(i,1);
if (t->getType() == TYPE_QUA){
u = -1 + 2*u;
v = -1 + 2*v;
@@ -350,8 +350,8 @@ double mesh_functional_distorsion_pN(MElement *t)
Ji(i) = mesh_functional_distorsion_2D(t,u,v);
}
- fullVector<double> Bi( jac->matrixLag2Bez.size1() );
- jac->matrixLag2Bez.mult(Ji,Bi);
+ fullVector<double> Bi( jac->getNumJacNodes() );
+ jac->lag2Bez(Ji,Bi);
return *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
}
@@ -421,17 +421,17 @@ double mesh_functional_distorsion_3D(MElement *t, double u, double v, double w)
double qmDistorsionOfMapping(MTetrahedron *t)
{
- const bezierBasis *jac = t->getJacobianFuncSpace()->bezier;
- fullVector<double>Ji(jac->points.size1());
- for (int i=0;i<jac->points.size1();i++){
- const double u = jac->points(i,0);
- const double v = jac->points(i,1);
- const double w = jac->points(i,2);
+ const JacobianBasis *jac = t->getJacobianFuncSpace();
+ fullVector<double>Ji(jac->getNumJacNodes());
+ for (int i=0;i<jac->getNumJacNodes();i++){
+ const double u = jac->getPoints()(i,0);
+ const double v = jac->getPoints()(i,1);
+ const double w = jac->getPoints()(i,2);
Ji(i) = mesh_functional_distorsion_3D(t,u,v,w);
}
- fullVector<double> Bi( jac->matrixLag2Bez.size1() );
- jac->matrixLag2Bez.mult(Ji,Bi);
+ fullVector<double> Bi( jac->getNumJacNodes() );
+ jac->lag2Bez(Ji,Bi);
/*
jac->matrixLag2Bez.print("Lag2Bez");
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 2a37328a38e4fa3e95fe220c64cb47fe0363db6e..a28fbe0554aa1dbb606efaa528461b2090caf931 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -935,27 +935,6 @@ static fullMatrix<double> generateSubDivisor
-static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &points, const nodalBasis *F)
-{
-
- fullMatrix<double> allDPsi;
- F->df(points, allDPsi);
-
- const int NBez = points.size1(), NLag = allDPsi.size1();
- jfs.gradShapeMatX.resize(NBez,NLag);
- jfs.gradShapeMatY.resize(NBez,NLag);
- jfs.gradShapeMatZ.resize(NBez,NLag);
- for (int i=0; i<NBez; i++)
- for (int j=0; j<NLag; j++) {
- jfs.gradShapeMatX(i,j) = allDPsi(j,3*i);
- jfs.gradShapeMatY(i,j) = allDPsi(j,3*i+1);
- jfs.gradShapeMatZ(i,j) = allDPsi(j,3*i+2);
- }
-
-}
-
-
-
std::map<int, bezierBasis> bezierBasis::_bbs;
const bezierBasis *bezierBasis::find(int tag)
{
@@ -967,6 +946,7 @@ const bezierBasis *bezierBasis::find(int tag)
B.order = MElement::OrderFromTag(tag);
if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
+ B.dim = 3;
B.numLagPts = 5;
B.points = gmshGeneratePointsPyramid(B.order,false);
B.matrixLag2Bez.resize(B.points.size1(),B.points.size1(),0.);
@@ -982,12 +962,14 @@ const bezierBasis *bezierBasis::find(int tag)
std::vector< fullMatrix<double> > subPoints;
switch (MElement::ParentTypeFromTag(tag)) {
case TYPE_PNT :
+ B.dim = 0;
B.numLagPts = 1;
B.exponents = generate1DExponents(0);
B.points = generate1DPoints(0);
dimSimplex = 0;
break;
case TYPE_LIN : {
+ B.dim = 1;
B.numLagPts = 2;
B.exponents = generate1DExponents(B.order);
B.points = generate1DPoints(B.order);
@@ -996,6 +978,7 @@ const bezierBasis *bezierBasis::find(int tag)
break;
}
case TYPE_TRI : {
+ B.dim = 2;
B.numLagPts = 3;
B.exponents = generateExponentsTriangle(B.order);
B.points = gmshGeneratePointsTriangle(B.order,false);
@@ -1004,6 +987,7 @@ const bezierBasis *bezierBasis::find(int tag)
break;
}
case TYPE_QUA : {
+ B.dim = 2;
B.numLagPts = 4;
B.exponents = generateExponentsQuad(B.order);
B.points = generatePointsQuad(B.order,false);
@@ -1013,6 +997,7 @@ const bezierBasis *bezierBasis::find(int tag)
break;
}
case TYPE_TET : {
+ B.dim = 3;
B.numLagPts = 4;
B.exponents = generateExponentsTetrahedron(B.order);
B.points = gmshGeneratePointsTetrahedron(B.order,false);
@@ -1021,6 +1006,7 @@ const bezierBasis *bezierBasis::find(int tag)
break;
}
case TYPE_PRI : {
+ B.dim = 3;
B.numLagPts = 6;
B.exponents = generateExponentsPrism(B.order);
B.points = generatePointsPrism(B.order, false);
@@ -1029,6 +1015,7 @@ const bezierBasis *bezierBasis::find(int tag)
break;
}
case TYPE_HEX : {
+ B.dim = 3;
B.numLagPts = 8;
B.exponents = generateExponentsHex(B.order);
B.points = generatePointsHex(B.order, false);
@@ -1038,6 +1025,7 @@ const bezierBasis *bezierBasis::find(int tag)
}
default : {
Msg::Error("Unknown function space %d: reverting to TET_1", tag);
+ B.dim = 3;
B.numLagPts = 4;
B.exponents = generateExponentsTetrahedron(0);
B.points = gmshGeneratePointsTetrahedron(0, false);
@@ -1057,30 +1045,18 @@ const bezierBasis *bezierBasis::find(int tag)
-static void generatePrimJac2Jac(JacobianBasis &jfs, const nodalBasis *primJacBasis)
+JacobianBasis::JacobianBasis(int tag)
{
- const int numJacNodes = jfs.bezier->points.size1(), numPrimJacSF = primJacBasis->getNumShapeFunctions();
- jfs.primJac2Jac.resize(numJacNodes,numPrimJacSF);
- primJacBasis->f(jfs.bezier->points,jfs.primJac2Jac);
-}
-
-
-std::map<int, JacobianBasis> JacobianBasis::_fs;
-const JacobianBasis *JacobianBasis::find(int tag)
-{
- std::map<int, JacobianBasis>::const_iterator it = _fs.find(tag);
- if (it != _fs.end()) return &it->second;
- JacobianBasis &J = _fs[tag];
int jacType, primJacType;
+
if (MElement::ParentTypeFromTag(tag) == TYPE_PYR) {
switch (tag) {
case MSH_PYR_5 : jacType = MSH_PYR_14; primJacType = MSH_PYR_14; break; // TODO: Order 1, Jac. "order" 2, check this
case MSH_PYR_14 : jacType = MSH_PYR_91; primJacType = MSH_PYR_14; break; // TODO: Order 2, Jac. "order" 5, check this
case MSH_PYR_30 : jacType = MSH_PYR_285; primJacType = MSH_PYR_14; break; // TODO: Order 3, Jac. "order" 8, check this
default :
- Msg::Error("Unknown Jacobian function space for element type %d: reverting to PYR_5", tag);
- J.bezier = bezierBasis::find(MSH_PYR_14);
+ Msg::Error("Unknown Jacobian function space for element type %d", tag);
break;
}
}
@@ -1103,54 +1079,74 @@ const JacobianBasis *JacobianBasis::find(int tag)
jacType = polynomialBasis::getTag(parentType, jacobianOrder, false);
primJacType = polynomialBasis::getTag(parentType, primJacobianOrder, false);
}
- J.bezier = bezierBasis::find(jacType);
+
+ // Store Bezier basis
+ bezier = bezierBasis::find(jacType);
+
+ // Store shape function gradients of mapping at Jacobian nodes
const nodalBasis *mapBasis = BasisFactory::create(tag);
- generateGradShapes(J, J.bezier->points, mapBasis);
+ fullMatrix<double> allDPsi;
+ mapBasis->df(getPoints(), allDPsi);
+ const int numJacNodes = getPoints().size1(), numMapNodes = allDPsi.size1();
+ gradShapeMatX.resize(numJacNodes,numMapNodes);
+ gradShapeMatY.resize(numJacNodes,numMapNodes);
+ gradShapeMatZ.resize(numJacNodes,numMapNodes);
+ for (int i=0; i<numJacNodes; i++)
+ for (int j=0; j<numMapNodes; j++) {
+ gradShapeMatX(i,j) = allDPsi(j,3*i);
+ gradShapeMatY(i,j) = allDPsi(j,3*i+1);
+ gradShapeMatZ(i,j) = allDPsi(j,3*i+2);
+ }
+
+ // Compute matrix for lifting from primary Jacobian basis to Jacobian basis
const nodalBasis *primJacBasis = BasisFactory::create(primJacType);
- generatePrimJac2Jac(J, primJacBasis);
- return &J;
+ const int numPrimJacSF = primJacBasis->getNumShapeFunctions();
+ matrixPrimJac2Jac.resize(numJacNodes,numPrimJacSF);
+ primJacBasis->f(getPoints(),matrixPrimJac2Jac);
+
}
-void JacobianBasis::getSignedJacobian(MElement *el, fullVector<double> &jacobian) const
+std::map<int, JacobianBasis*> JacobianBasis::_fs;
+const JacobianBasis *JacobianBasis::find(int tag)
{
- const int numVertices = gradShapeMatX.size2(), numJacNodes = gradShapeMatX.size1();
+ std::map<int, JacobianBasis*>::const_iterator it = _fs.find(tag);
+ if (it != _fs.end()) return it->second;
+
+ JacobianBasis *B = new JacobianBasis(tag);
+ _fs.insert(std::make_pair(tag, B));
+ return B;
+
+}
+
+
- switch (el->getDim()) {
+void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
+{
+
+ switch (bezier->dim) {
case 1 : {
+ const int numVertices = gradShapeMatX.size2();
fullVector<double> nodesX(numVertices);
- for (int i = 0; i < numVertices; i++) {
- nodesX(i) = el->getShapeFunctionNode(i)->x();
- }
+ nodesXYZ.copyOneColumn(nodesX,0);
gradShapeMatX.mult(nodesX, jacobian);
break;
}
case 2 : {
- fullMatrix<double> nodesXY(numVertices,2);
- for (int i = 0; i < numVertices; i++) {
- MVertex *v = el->getShapeFunctionNode(i);
- nodesXY(i,0) = v->x();
- nodesXY(i,1) = v->y();
- }
+ const int numJacNodes = gradShapeMatX.size1();
fullMatrix<double> dxydX(numJacNodes,2), dxydY(numJacNodes,2);
- gradShapeMatX.mult(nodesXY, dxydX);
- gradShapeMatY.mult(nodesXY, dxydY);
+ gradShapeMatX.mult(nodesXYZ, dxydX);
+ gradShapeMatY.mult(nodesXYZ, dxydY);
for (int i = 0; i < numJacNodes; i++) jacobian(i) = dxydX(i,0)*dxydY(i,1)-dxydX(i,1)*dxydY(i,0);
break;
}
case 3 : {
- fullMatrix<double> nodesXYZ(numVertices,3);
- for (int i = 0; i < numVertices; i++) {
- MVertex *v = el->getShapeFunctionNode(i);
- nodesXYZ(i,0) = v->x();
- nodesXYZ(i,1) = v->y();
- nodesXYZ(i,2) = v->z();
- }
+ const int numJacNodes = gradShapeMatX.size1();
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
gradShapeMatY.mult(nodesXYZ, dxyzdY);
@@ -1168,3 +1164,51 @@ void JacobianBasis::getSignedJacobian(MElement *el, fullVector<double> &jacobian
}
}
+
+
+
+void JacobianBasis::getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
+{
+
+ switch (bezier->dim) {
+
+ case 1 : {
+ gradShapeMatX.mult(nodesX, jacobian);
+ break;
+ }
+
+ case 2 : {
+ const int numJacNodes = gradShapeMatX.size1(), numEl = nodesX.size2();
+ fullMatrix<double> dxdX(numJacNodes,numEl), dydX(numJacNodes,numEl);
+ fullMatrix<double> dxdY(numJacNodes,numEl), dydY(numJacNodes,numEl);
+ gradShapeMatX.mult(nodesX, dxdX); gradShapeMatX.mult(nodesY, dydX);
+ gradShapeMatY.mult(nodesX, dxdY); gradShapeMatY.mult(nodesY, dydY);
+ for (int iEl = 0; iEl < numEl; iEl++)
+ for (int i = 0; i < numJacNodes; i++)
+ jacobian(i,iEl) = dxdX(i,iEl)*dydY(i,iEl)-dydX(i,iEl)*dxdY(i,iEl);
+ break;
+ }
+
+ case 3 : {
+ const int numJacNodes = gradShapeMatX.size1(), numEl = nodesX.size2();
+ fullMatrix<double> dxdX(numJacNodes,numEl), dydX(numJacNodes,numEl), dzdX(numJacNodes,numEl);
+ fullMatrix<double> dxdY(numJacNodes,numEl), dydY(numJacNodes,numEl), dzdY(numJacNodes,numEl);
+ fullMatrix<double> dxdZ(numJacNodes,numEl), dydZ(numJacNodes,numEl), dzdZ(numJacNodes,numEl);
+ gradShapeMatX.mult(nodesX, dxdX); gradShapeMatX.mult(nodesY, dydX); gradShapeMatX.mult(nodesZ, dzdX);
+ gradShapeMatY.mult(nodesX, dxdY); gradShapeMatY.mult(nodesY, dydY); gradShapeMatY.mult(nodesZ, dzdY);
+ gradShapeMatZ.mult(nodesX, dxdZ); gradShapeMatZ.mult(nodesY, dydZ); gradShapeMatZ.mult(nodesZ, dzdZ);
+ for (int iEl = 0; iEl < numEl; iEl++)
+ for (int i = 0; i < numJacNodes; i++)
+ jacobian(i,iEl) = dxdX(i,iEl)*dydY(i,iEl)*dzdZ(i,iEl)
+ + dxdY(i,iEl)*dydZ(i,iEl)*dzdX(i,iEl)
+ + dydX(i,iEl)*dzdY(i,iEl)*dxdZ(i,iEl)
+ - dxdZ(i,iEl)*dydY(i,iEl)*dzdX(i,iEl)
+ - dxdY(i,iEl)*dydX(i,iEl)*dzdZ(i,iEl)
+ - dydZ(i,iEl)*dzdY(i,iEl)*dxdX(i,iEl);
+ break;
+ }
+
+ }
+
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index c7cf23c4785deb236ece3673edfdffac2937f93a..d5761cdde3410246a396dcabc37d24a5c8fc1894 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -16,7 +16,7 @@ class bezierBasis {
private :
static std::map<int, bezierBasis> _bbs;
public :
- int order;
+ int dim, order;
int numLagPts;
int numDivisions;
// the 'numLagPts' first exponents are related to 'Lagrangian' values
@@ -29,13 +29,34 @@ class bezierBasis {
class JacobianBasis {
private:
- static std::map<int, JacobianBasis> _fs;
- public :
+ static std::map<int, JacobianBasis*> _fs;
const bezierBasis *bezier;
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
- fullMatrix<double> primJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
+ fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
+ public :
static const JacobianBasis *find(int);
- void getSignedJacobian(MElement *el, fullVector<double> &jacobian) const;
+ JacobianBasis(int tag);
+ inline int getNumJacNodes() const { return gradShapeMatX.size1(); }
+ inline int getNumMapNodes() const { return gradShapeMatX.size2(); }
+ inline const fullMatrix<double> &getPoints() const { return bezier->points; }
+ inline int getNumDivisions() const { return bezier->numDivisions; }
+ inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
+ inline int getNumLagPts() const { return bezier->numLagPts; }
+ void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
+ void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
+ inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
+ bezier->matrixLag2Bez.mult(jac,bez);
+ }
+ inline void lag2Bez(const fullMatrix<double> &jac, fullMatrix<double> &bez) const {
+ bezier->matrixLag2Bez.mult(jac,bez);
+ }
+ inline void primJac2Jac(const fullVector<double> &primJac, fullVector<double> &jac) const {
+ matrixPrimJac2Jac.mult(primJac,jac);
+ }
+ inline void subDivisor(const fullVector<double> &bez, fullVector<double> &result) const {
+ bezier->subDivisor.mult(bez,result);
+ }
};
#endif
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index 87036b090cb269966054012fee95fef3cdb3ded4..efbc446b1a5ed5457f7fe53191cd70eae42d60fd 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -228,23 +228,39 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
- const bezierBasis *bb = jfs->bezier;
-
- int numSamplingPt = bb->points.size1();
+ const int numSamplingPt = jfs->getNumJacNodes(), numSamplingPt1 = jfs1->getNumJacNodes();
+ const int numMapNodes = jfs->getNumMapNodes(), numMapNodes1 = jfs1->getNumMapNodes();
+ const int dim = el[0]->getDim();
#ifdef _ANALYSECURVEDMESH_BLAS_
fullMatrix<double> jacobianB(numSamplingPt, numEl);
fullMatrix<double> jacBezB(numSamplingPt, numEl);
- fullVector<double> jac1B(jfs1->bezier->points.size1(), numEl);
+ fullMatrix<double> jac1B(numSamplingPt1, numEl);
fullVector<double> jacBez, jacobian, jac1;
- jfs->getSignedJacobian(el, jacobianB);
- jfs1->getSignedJacobian(el, jac1B);
- bb->matrixLag2Bez.mult(jacobianB, jacBezB);
+ fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
+ fullMatrix<double> nodesX1(numMapNodes, numEl), nodesY1(numMapNodes, numEl), nodesZ1(numMapNodes, numEl);
+ for (int k = 0; k < numEl; ++k)
+ for (int i = 0; i < numMapNodes; ++i)
+ {
+ MVertex *v = (el[k])->getShapeFunctionNode(i);
+ nodesX(i,k) = v->x();
+ nodesY(i,k) = v->y();
+ nodesZ(i,k) = v->z();
+ if (i < numMapNodes1) {
+ nodesX1(i,k) = nodesX(i,k);
+ nodesY1(i,k) = nodesY(i,k);
+ nodesZ1(i,k) = nodesZ(i,k);
+ }
+ }
+
+ jfs->getSignedJacobian(nodesX, nodesY, nodesZ, jacobianB);
+ jfs1->getSignedJacobian(nodesX1, nodesY1, nodesZ1, jac1B);
+ jfs->lag2Bez(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
- fullVector<double> jac1(jfs1->bezier->points.size1());
+ fullVector<double> jac1(numSamplingPt1);
#endif
for (int k = 0; k < numEl; ++k) {
@@ -254,8 +270,11 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
jacobian.setAsProxy(jacobianB, k);
jac1.setAsProxy(jac1B, k);
#else
- jfs->getSignedJacobian(el[k], jacobian);
- jfs1->getSignedJacobian(el[k], jac1);
+ fullMatrix<double> nodesXYZ(numMapNodes,dim), nodesXYZ1(numMapNodes1,dim);
+ el[k]->getNodesCoord(nodesXYZ);
+ nodesXYZ1.copy(nodesXYZ,0,numMapNodes1,0,dim,0,0);
+ jfs->getSignedJacobian(nodesXYZ,jacobian);
+ jfs1->getSignedJacobian(nodesXYZ1,jac1);
#endif
// AmJ : avgJ is not the average Jac for quad, prism or hex
@@ -281,7 +300,7 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
}
#ifndef _ANALYSECURVEDMESH_BLAS_
- bb->matrixLag2Bez.mult(jacobian, jacBez);
+ jfs->lag2Bez(jacobian, jacBez);
#endif
for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak * avgJ; ++i);
@@ -296,7 +315,7 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
continue;
}
else {
- int result = subDivision(bb, jacBez, _maxDepth-1);
+ int result = subDivision(jfs, jacBez, _maxDepth-1);
if (result < 0) {
invalids.push_back(el[k]);
++_numInvalid;
@@ -313,16 +332,16 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
}
}
-int GMSH_AnalyseCurvedMeshPlugin::subDivision(const bezierBasis *bb,
+int GMSH_AnalyseCurvedMeshPlugin::subDivision(const JacobianBasis *jb,
const fullVector<double> &jacobian,
int depth)
{
- fullVector<double> newJacobian(bb->subDivisor.size1());
- bb->subDivisor.mult(jacobian, newJacobian);
+ fullVector<double> newJacobian(jb->getNumSubNodes());
+ jb->subDivisor(jacobian, newJacobian);
- for (int i = 0; i < bb->numDivisions; i++)
- for (int j = 0; j < bb->numLagPts; j++)
- if (newJacobian(i * bb->points.size1() + j) <= _jacBreak)
+ for (int i = 0; i < jb->getNumDivisions(); i++)
+ for (int j = 0; j < jb->getNumLagPts(); j++)
+ if (newJacobian(i * jb->getNumJacNodes() + j) <= _jacBreak)
return -1;
int i = 0;
@@ -338,9 +357,9 @@ int GMSH_AnalyseCurvedMeshPlugin::subDivision(const bezierBasis *bb,
std::vector<int> negTag, posTag;
bool zeroTag = false;
- for (int i = 0; i < bb->numDivisions; i++) {
+ for (int i = 0; i < jb->getNumDivisions(); i++) {
subJacobian.setAsProxy(newJacobian, i * jacobian.size(), jacobian.size());
- int tag = subDivision(bb, subJacobian, depth-1);
+ int tag = subDivision(jb, subJacobian, depth-1);
if (tag < 0)
negTag.push_back(tag);
@@ -433,25 +452,42 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
- const bezierBasis *bb = jfs->bezier;
- int numSamplingPt = bb->points.size1();
+ const int numSamplingPt = jfs->getNumJacNodes(), numSamplingPt1 = jfs1->getNumJacNodes();
+ const int numMapNodes = jfs->getNumMapNodes(), numMapNodes1 = jfs1->getNumMapNodes();
+ const int dim = el[0]->getDim();
#ifdef _ANALYSECURVEDMESH_BLAS_
fullMatrix<double> jacobianB(numSamplingPt, numEl);
fullMatrix<double> jacBezB(numSamplingPt, numEl);
- fullVector<double> jac1B(jfs1->bezier->points.size1(), numEl);
+ fullMatrix<double> jac1B(numSamplingPt1, numEl);
fullVector<double> jacBez, jacobian, jac1;
- jfs->getSignedJacobian(el, jacobianB);
- jfs1->getSignedJacobian(el, jac1B);
- bb->matrixLag2Bez.mult(jacobianB, jacBezB);
+ fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
+ fullMatrix<double> nodesX1(numMapNodes1, numEl), nodesY1(numMapNodes1, numEl), nodesZ1(numMapNodes1, numEl);
+ for (int k = 0; k < numEl; ++k)
+ for (int i = 0; i < numMapNodes; ++i)
+ {
+ MVertex *v = (el[k])->getShapeFunctionNode(i);
+ nodesX(i,k) = v->x();
+ nodesY(i,k) = v->y();
+ nodesZ(i,k) = v->z();
+ if (i < numMapNodes1) {
+ nodesX1(i,k) = nodesX(i,k);
+ nodesY1(i,k) = nodesY(i,k);
+ nodesZ1(i,k) = nodesZ(i,k);
+ }
+ }
+
+ jfs->getSignedJacobian(nodesX, nodesY, nodesZ, jacobianB);
+ jfs1->getSignedJacobian(nodesX, nodesY, nodesZ, jac1B);
+ jfs->lag2Bez(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
- fullVector<double> jac1(jfs1->bezier->points.size1());
+ fullVector<double> jac1(numSamplingPt1);
#endif
- fullVector<double> subJacBez(bb->subDivisor.size1());
+ fullVector<double> subJacBez(jfs->getNumSubNodes());
_min_Javg = 1.7e308;
_max_Javg = -1.7e308;
@@ -468,9 +504,12 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
jacobian.setAsProxy(jacobianB, k);
jac1.setAsProxy(jac1B, k);
#else
- jfs->getSignedJacobian(el[k], jacobian);
- jfs1->getSignedJacobian(el[k], jac1);
- bb->matrixLag2Bez.mult(jacobian, jacBez);
+ fullMatrix<double> nodesXYZ(numMapNodes,dim), nodesXYZ1(numMapNodes1,dim);
+ el[k]->getNodesCoord(nodesXYZ);
+ nodesXYZ1.copy(nodesXYZ,0,numMapNodes1,0,dim,0,0);
+ jfs->getSignedJacobian(nodesXYZ,jacobian);
+ jfs1->getSignedJacobian(nodesXYZ1,jac1);
+ jfs->lag2Bez(jacobian, jacBez);
#endif
// AmJ : avgJ is not the average Jac for quad, prism or hex
@@ -522,7 +561,7 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
pqMin.pop();
delete bj;
- for (int i = 0; i < bb->numDivisions; i++) {
+ for (int i = 0; i < jfs->getNumDivisions(); i++) {
jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
bj = new BezierJacobian(jacBez, jfs, currentDepth);
pqMin.push(bj);
@@ -555,7 +594,7 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
pqMax.pop();
delete bj;
- for (int i = 0; i < bb->numDivisions; i++) {
+ for (int i = 0; i < jfs->getNumDivisions(); i++) {
jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
bj = new BezierJacobian(jacBez, jfs, currentDepth);
pqMax.push(bj);
@@ -690,7 +729,7 @@ BezierJacobian::BezierJacobian(fullVector<double> &v, const JacobianBasis *jfs,
_minJ = _maxJ = v(0);
int i = 1;
- for (; i < jfs->bezier->numLagPts; i++) {
+ for (; i < jfs->getNumLagPts(); i++) {
if (_minJ > v(i)) _minJ = v(i);
if (_maxJ < v(i)) _maxJ = v(i);
}
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index fbbb2c76a9152c2dfd780f9f9617c3c450ce13c8..73b933f2344d22da7e896b2a94addb73172a93a3 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -49,7 +49,7 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
void checkValidity(int toDo);
void computeMinMax(std::map<int, std::vector<double> > *data = 0);
void computeMinMax(MElement *const *, int numEl, std::map<int, std::vector<double> > *data = 0);
- int subDivision(const bezierBasis*, const fullVector<double>&, int depth);
+ int subDivision(const JacobianBasis *, const fullVector<double>&, int depth);
void hideValid_ShowInvalid(std::vector<MElement*> &invalids);
};
@@ -72,7 +72,7 @@ private:
public:
BezierJacobian(fullVector<double> &, const JacobianBasis *, int depth);
void subDivisions(fullVector<double> &vect) const
- {_jfs->bezier->subDivisor.mult(_jacBez, vect);}
+ {_jfs->subDivisor(_jacBez, vect);}
inline int depth() const {return _depthSub;}
inline double minJ() const {return _minJ;}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index a74db6c0ea3594c8898f6f744d071c17404cfd8e..12c2c6dedfcacdf4ee4e2644021c47319bf8c1fb 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -9,14 +9,13 @@
std::map<int, fullMatrix<double> > Mesh::_gradShapeFunctions;
-std::map<int, fullMatrix<double> > Mesh::_lag2Bez;
-fullMatrix<double> Mesh::computeGSF(const nodalBasis *lagrange, const bezierBasis *bezier)
+fullMatrix<double> Mesh::computeGSF(const nodalBasis *lagrange, const JacobianBasis *jac)
{
// bezier points are defined in the [0,1] x [0,1] quad
- fullMatrix<double> bezierPoints = bezier->points;
+ fullMatrix<double> bezierPoints = jac->getPoints();
if (lagrange->parentType == TYPE_QUA) {
for (int i = 0; i < bezierPoints.size1(); ++i) {
bezierPoints(i, 0) = -1 + 2 * bezierPoints(i, 0);
@@ -71,12 +70,10 @@ Mesh::Mesh(GEntity *ge, const std::set<MElement*> &els, std::set<MVertex*> &toFi
MElement *el = *it;
_el[iEl] = el;
const nodalBasis *lagrange = el->getFunctionSpace();
- const bezierBasis *bezier = JacobianBasis::find(lagrange->type)->bezier;
- if (_lag2Bez.find(lagrange->type) == _lag2Bez.end()) {
- _gradShapeFunctions[lagrange->type] = computeGSF(lagrange, bezier);
- _lag2Bez[lagrange->type] = bezier->matrixLag2Bez;
- }
- _nBezEl[iEl] = bezier->points.size1();
+ const JacobianBasis *jac = JacobianBasis::find(lagrange->type);
+ if (_gradShapeFunctions.find(lagrange->type) == _gradShapeFunctions.end())
+ _gradShapeFunctions[lagrange->type] = computeGSF(lagrange, jac);
+ _nBezEl[iEl] = jac->getNumJacNodes();
_nNodEl[iEl] = lagrange->points.size1();
for (int iVEl = 0; iVEl < lagrange->points.size1(); iVEl++) {
MVertex *vert = el->getVertex(iVEl);
@@ -375,8 +372,8 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<
// gradScaledJac(iEl,OLD);
// scaledJac(iEl,OLD);
+ const JacobianBasis *jacBasis = JacobianBasis::find(_el[iEl]->getTypeForMSH());
fullMatrix<double> &gsf = _gradShapeFunctions[_el[iEl]->getTypeForMSH()];
- const fullMatrix<double> &l2b = _lag2Bez[_el[iEl]->getTypeForMSH()];
const int nbBez = _nBezEl[iEl];
const int nbNod = _nNodEl[iEl];
fullMatrix<double> JDJ (nbBez,3*nbNod+1);
@@ -428,7 +425,7 @@ void Mesh::scaledJacAndGradients(int iEl, std::vector<double> &sJ , std::vector<
}
// (N_b x N_b) x (N_b x 3*N_n + 1)
- l2b.mult(JDJ,BDB);
+ jacBasis->lag2Bez(JDJ,BDB);
// the scaled jacobian
// printf("ELEMENT %d\n",iEl);
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.h b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
index 87bdd18ffd7ce7cc4fb4e9d97ec6e9aab178eaf2..a7fa37a76848e2fb834036ea05b8d957423f3c67 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.h
@@ -80,12 +80,11 @@ private:
bool projJac; // Using "projected" Jacobians or not
static std::map<int, fullMatrix<double> > _gradShapeFunctions; // gradients of shape functions at Bezier points
- static std::map<int, fullMatrix<double> > _lag2Bez; // gradients of shape functions at Bezier points
int addVert(MVertex* vert);
int addFreeVert(MVertex* vert, const int iV, const int nPCV, std::set<MVertex*> &toFix);
SVector3 getNormalEl(int iEl);
- static fullMatrix<double> computeGSF(const nodalBasis *lagrange, const bezierBasis *bezier);
+ static fullMatrix<double> computeGSF(const nodalBasis *lagrange, const JacobianBasis *jac);
static inline int indJB2DBase(int nNod, int l, int i, int j) { return (l*nNod+i)*nNod+j; }
inline int indJB2D(int iEl, int l, int i, int j) { return indJB2DBase(_nNodEl[iEl],l,i,j); }
static inline int indJB3DBase(int nNod, int l, int i, int j, int m) { return ((l*nNod+i)*nNod+j)*nNod+m; }