diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 5f0d00f67fba8fcbfe5efa42e42f7414ebcd4a67..7a89cdff8c44fa5cf6574a4bb57640f3d1db82c2 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -112,28 +112,22 @@ void MElement::scaledJacRange(double &jmin, double &jmax)
{
jmin = jmax = 1.0;
#if defined(HAVE_MESH)
- extern double mesh_functional_distorsion_2D(MElement*,double,double);
- extern double mesh_functional_distorsion_3D(MElement*,double,double,double);
if (getPolynomialOrder() == 1) return;
- const bezierBasis *jac = 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);
- if (getDim() == 2) {
- if (getType() == TYPE_QUA) {
- u = -1 + 2*u;
- v = -1 + 2*v;
- }
- Ji(i) = mesh_functional_distorsion_2D(this,u,v);
- }
- else {
- double w = jac->points(i,2);
- Ji(i) = mesh_functional_distorsion_3D(this,u,v,w);
- }
- }
- fullVector<double> Bi( jac->matrixLag2Bez.size1() );
- jac->matrixLag2Bez.mult(Ji,Bi);
+ 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
+ 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);
jmin = *std::min_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
jmax = *std::max_element(Bi.getDataPtr(),Bi.getDataPtr()+Bi.size());
#endif
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 9029b64772a471e223c2cb82f842d81f1b9cf73c..c6f3a61af706a66e862f80d6dfb5f43d0ba5673b 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -182,7 +182,7 @@ class MElement
virtual double rhoShapeMeasure();
virtual double gammaShapeMeasure(){ return 0.; }
virtual double etaShapeMeasure(){ return 0.; }
- virtual double distoShapeMeasure(){ return 1.0; }
+ virtual double distoShapeMeasure(){ double jmin,jmax; scaledJacRange(jmin,jmax); return jmin; }
virtual double angleShapeMeasure() { return 1.0; }
virtual void scaledJacRange(double &jmin, double &jmax);
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index f4e42ef6ca7a4fb60502618e5476fce80723cfd7..b8b1ab8373ed42df18bef44668f9a2832885c7d0 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -82,16 +82,6 @@ const JacobianBasis* MPyramid::getJacobianFuncSpace(int o) const
MPyramidN::~MPyramidN() {}
-double MPyramidN::distoShapeMeasure()
-{
-#if defined(HAVE_MESH)
- _disto = 0.0; //qmDistorsionOfMapping(this);
-#else
- _disto = 0.;
-#endif
- return _disto;
-}
-
int MPyramidN::getNumEdgesRep(){ return 8 * CTX::instance()->mesh.numSubEdges; }
void MPyramidN::getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 836bff3b3595ff38e601f70904edd91f7e8efeb0..ffae41f7d1cfd15a0983c2339711bb155be5446d 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -436,7 +436,6 @@ class MPyramidN : public MPyramid {
~MPyramidN();
- virtual double distoShapeMeasure();
virtual int getPolynomialOrder() const { return _order; }
virtual int getNumVertices() const { return 5 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 5 ? _v[num] : _vs[num - 5]; }
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 1dde8a9e9a5cd385dfd67fde9dc7241e74948d00..6265eb91fe33eab526a17eaaf100ec0350d1803b 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -119,16 +119,6 @@ double MTetrahedron::distoShapeMeasure()
#endif
}
-double MTetrahedronN::distoShapeMeasure()
-{
-#if defined(HAVE_MESH)
- _disto = qmDistorsionOfMapping(this);
-#else
- _disto = 0.;
-#endif
- return _disto;
-}
-
double MTetrahedron::getInnerRadius()
{
// radius of inscribed sphere = 3 * Volume / sum(Area_i)
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index 2df127b1640bc3c12d7ebd9bae6fd3b709705576..d0e5e9cae82f934bf64b3baa9be8ccbfd10d7066 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -213,7 +213,6 @@ class MTetrahedron10 : public MTetrahedron {
for(int i = 0; i < 6; i++) _vs[i]->setPolynomialOrder(2);
}
~MTetrahedron10(){}
- virtual double distoShapeMeasure() {double jmin,jmax; scaledJacRange(jmin,jmax); return jmin;}
virtual void scaledJacRange(double &jmin, double &jmax);
virtual int getPolynomialOrder() const { return 2; }
virtual int getNumVertices() const { return 10; }
@@ -315,22 +314,20 @@ class MTetrahedronN : public MTetrahedron {
protected:
std::vector<MVertex *> _vs;
const char _order;
- double _disto;
public:
MTetrahedronN(MVertex *v0, MVertex *v1, MVertex *v2, MVertex *v3,
const std::vector<MVertex*> &v, char order, int num=0, int part=0)
- : MTetrahedron(v0, v1, v2, v3, num, part) , _vs (v), _order(order),_disto(-1.e22)
+ : MTetrahedron(v0, v1, v2, v3, num, part) , _vs (v), _order(order)
{
for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
}
MTetrahedronN(const std::vector<MVertex*> &v, char order, int num=0, int part=0)
- : MTetrahedron(v[0], v[1], v[2], v[3], num, part) , _order(order),_disto(-1.e22)
+ : MTetrahedron(v[0], v[1], v[2], v[3], num, part) , _order(order)
{
for(unsigned int i = 4; i < v.size(); i++) _vs.push_back(v[i]);
for(unsigned int i = 0; i < _vs.size(); i++) _vs[i]->setPolynomialOrder(_order);
}
~MTetrahedronN(){}
- virtual double distoShapeMeasure();
virtual int getPolynomialOrder() const { return _order; }
virtual int getNumVertices() const { return 4 + _vs.size(); }
virtual MVertex *getVertex(int num){ return num < 4 ? _v[num] : _vs[num - 4]; }
diff --git a/Numeric/BergotBasis.h b/Numeric/BergotBasis.h
index e62e43828269978a4c19217870fbafc74294e2c0..95dacaf44036590871acc8a43a32e4c2d3c1a647 100644
--- a/Numeric/BergotBasis.h
+++ b/Numeric/BergotBasis.h
@@ -11,6 +11,11 @@
#include "jacobiPolynomials.h"
#include "legendrePolynomials.h"
+
+// Basis functios for pyramidal elements
+// cf. M. Bergot, G. Cohen, M. Durufle, HIGHER-ORDER FINITE ELEMENTS FOR HYBRID MESHES USING NEW
+// NODAL PYRAMIDAL ELEMENTS, J. Sci. Comput. 42, 3 (2010) 345-381", DOI: 10.1007/s10915-009-9334-9
+
class BergotBasis {
public:
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index ab0b938ec0bdf6a8d56668d65f325fa823a3f852..2a37328a38e4fa3e95fe220c64cb47fe0363db6e 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -11,6 +11,7 @@
#include "pyramidalBasis.h"
#include "pointsGenerators.h"
#include "BasisFactory.h"
+#include "MElement.h"
// Bezier Exponents
static fullMatrix<double> generate1DExponents(int order)
@@ -934,98 +935,7 @@ static fullMatrix<double> generateSubDivisor
-static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &points,
- const fullMatrix<double> &monomials, const fullMatrix<double> &coefficients)
-{
-
- int nbPts = points.size1();
- int nbDofs = monomials.size1();
- int dim = points.size2();
-
- switch (dim) {
- case 3 :
- jfs.gradShapeMatZ.resize(nbPts, nbDofs, true);
- case 2 :
- jfs.gradShapeMatY.resize(nbPts, nbDofs, true);
- case 1 :
- jfs.gradShapeMatX.resize(nbPts, nbDofs, true);
- break;
- default :
- return;
- }
-
- double dx, dy, dz;
-
- switch (dim) {
- case 3 :
- for (int i = 0; i < nbDofs; i++) {
- for (int k = 0; k < nbPts; k++) {
-
- if ((int) monomials(i, 0) > 0) {
- dx = pow( points(k, 0), monomials(i, 0)-1 ) * monomials(i, 0)
- * pow( points(k, 1), monomials(i, 1) )
- * pow( points(k, 2), monomials(i, 2) );
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
- }
- if ((int) monomials(i, 1) > 0.) {
- dy = pow( points(k, 0), monomials(i, 0) )
- * pow( points(k, 1), monomials(i, 1)-1 ) * monomials(i, 1)
- * pow( points(k, 2), monomials(i, 2) );
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatY(k, j) += coefficients(j, i) * dy;
- }
- if ((int) monomials(i, 2) > 0.) {
- dz = pow( points(k, 0), monomials(i, 0) )
- * pow( points(k, 1), monomials(i, 1) )
- * pow( points(k, 2), monomials(i, 2)-1 ) * monomials(i, 2);
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatZ(k, j) += coefficients(j, i) * dz;
- }
- }
- }
- return;
-
- case 2 :
- for (int i = 0; i < nbDofs; i++) {
- for (int k = 0; k < nbPts; k++) {
-
- if ((int) monomials(i, 0) > 0) {
- dx = pow( points(k, 0), (int) monomials(i, 0)-1 ) * monomials(i, 0)
- * pow( points(k, 1), (int) monomials(i, 1) );
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
- }
- if ((int) monomials(i, 1) > 0) {
- dy = pow( points(k, 0), (int) monomials(i, 0) )
- * pow( points(k, 1), (int) monomials(i, 1)-1 ) * monomials(i, 1);
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatY(k, j) += coefficients(j, i) * dy;
- }
- }
- }
- return;
-
- case 1 :
- for (int i = 0; i < nbDofs; i++) {
- for (int k = 0; k < nbPts; k++) {
-
- if ((int) monomials(i, 0) > 0) {
- dx = pow( points(k, 0), (int) monomials(i, 0)-1 ) * monomials(i, 0);
- for (int j = 0; j < nbDofs; j++)
- jfs.gradShapeMatX(k, j) += coefficients(j, i) * dx;
- }
- }
- }
- break;
- }
- return;
-
-}
-
-
-
-static void generateGradShapesPyramid(JacobianBasis &jfs, const fullMatrix<double> &points, const pyramidalBasis *F)
+static void generateGradShapes(JacobianBasis &jfs, const fullMatrix<double> &points, const nodalBasis *F)
{
fullMatrix<double> allDPsi;
@@ -1045,6 +955,7 @@ static void generateGradShapesPyramid(JacobianBasis &jfs, const fullMatrix<doubl
}
+
std::map<int, bezierBasis> bezierBasis::_bbs;
const bezierBasis *bezierBasis::find(int tag)
{
@@ -1064,14 +975,12 @@ const bezierBasis *bezierBasis::find(int tag)
B.matrixLag2Bez(i,i) = 1.;
B.matrixBez2Lag(i,i) = 1.;
}
-// B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, F->order, dimSimplex);
-// B.numDivisions = (int) pow(2., (int) B.points.size2());
+ // TODO: Implement subdidivsor
}
else {
- const polynomialBasis *F = (polynomialBasis*)BasisFactory::create(tag);
int dimSimplex;
std::vector< fullMatrix<double> > subPoints;
- switch (F->parentType) {
+ switch (MElement::ParentTypeFromTag(tag)) {
case TYPE_PNT :
B.numLagPts = 1;
B.exponents = generate1DExponents(0);
@@ -1080,51 +989,51 @@ const bezierBasis *bezierBasis::find(int tag)
break;
case TYPE_LIN : {
B.numLagPts = 2;
- B.exponents = generate1DExponents(F->order);
- B.points = generate1DPoints(F->order);
+ B.exponents = generate1DExponents(B.order);
+ B.points = generate1DPoints(B.order);
dimSimplex = 0;
subPoints = generateSubPointsLine(0);
break;
}
case TYPE_TRI : {
B.numLagPts = 3;
- B.exponents = generateExponentsTriangle(F->order);
- B.points = gmshGeneratePointsTriangle(F->order,false);
+ B.exponents = generateExponentsTriangle(B.order);
+ B.points = gmshGeneratePointsTriangle(B.order,false);
dimSimplex = 2;
- subPoints = generateSubPointsTriangle(F->order);
+ subPoints = generateSubPointsTriangle(B.order);
break;
}
case TYPE_QUA : {
B.numLagPts = 4;
- B.exponents = generateExponentsQuad(F->order);
- B.points = generatePointsQuad(F->order,false);
+ B.exponents = generateExponentsQuad(B.order);
+ B.points = generatePointsQuad(B.order,false);
dimSimplex = 0;
- subPoints = generateSubPointsQuad(F->order);
+ subPoints = generateSubPointsQuad(B.order);
// B.points.print("points");
break;
}
case TYPE_TET : {
B.numLagPts = 4;
- B.exponents = generateExponentsTetrahedron(F->order);
- B.points = gmshGeneratePointsTetrahedron(F->order,false);
+ B.exponents = generateExponentsTetrahedron(B.order);
+ B.points = gmshGeneratePointsTetrahedron(B.order,false);
dimSimplex = 3;
- subPoints = generateSubPointsTetrahedron(F->order);
+ subPoints = generateSubPointsTetrahedron(B.order);
break;
}
case TYPE_PRI : {
B.numLagPts = 6;
- B.exponents = generateExponentsPrism(F->order);
- B.points = generatePointsPrism(F->order, false);
+ B.exponents = generateExponentsPrism(B.order);
+ B.points = generatePointsPrism(B.order, false);
dimSimplex = 2;
- subPoints = generateSubPointsPrism(F->order);
+ subPoints = generateSubPointsPrism(B.order);
break;
}
case TYPE_HEX : {
B.numLagPts = 8;
- B.exponents = generateExponentsHex(F->order);
- B.points = generatePointsHex(F->order, false);
+ B.exponents = generateExponentsHex(B.order);
+ B.points = generatePointsHex(B.order, false);
dimSimplex = 0;
- subPoints = generateSubPointsHex(F->order);
+ subPoints = generateSubPointsHex(B.order);
break;
}
default : {
@@ -1137,53 +1046,125 @@ const bezierBasis *bezierBasis::find(int tag)
break;
}
}
- B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.points, F->order, dimSimplex);
+ B.matrixBez2Lag = generateBez2LagMatrix(B.exponents, B.points, B.order, dimSimplex);
B.matrixBez2Lag.invert(B.matrixLag2Bez);
- B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, F->order, dimSimplex);
+ B.subDivisor = generateSubDivisor(B.exponents, subPoints, B.matrixLag2Bez, B.order, dimSimplex);
B.numDivisions = (int) pow(2., (int) B.points.size2());
}
return &B;
}
+
+
+static void generatePrimJac2Jac(JacobianBasis &jfs, const nodalBasis *primJacBasis)
+{
+ 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 : J.bezier = bezierBasis::find(MSH_PYR_14); break; // TODO: Order 1, Jac. "order" 2, check this
- case MSH_PYR_14 : J.bezier = bezierBasis::find(MSH_PYR_91); break; // TODO: Order 2, Jac. "order" 5, check this
- case MSH_PYR_30 : J.bezier = bezierBasis::find(MSH_PYR_285); break; // TODO: Order 3, Jac. "order" 8, check this
+ 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);
break;
}
- pyramidalBasis *F = (pyramidalBasis*)BasisFactory::create(tag);
- generateGradShapesPyramid(J, J.bezier->points, F);
}
else {
const int parentType = MElement::ParentTypeFromTag(tag), order = MElement::OrderFromTag(tag);
- int jacobianOrder;
+ int jacobianOrder, primJacobianOrder;
switch (parentType) {
- case TYPE_PNT : jacobianOrder = 0; break;
- case TYPE_LIN : jacobianOrder = order - 1; break;
- case TYPE_TRI : jacobianOrder = 2 * (order - 1); break;
- case TYPE_QUA : jacobianOrder = 2 * order - 1; break;
- case TYPE_TET : jacobianOrder = 3 * (order - 1); break;
- case TYPE_PRI : jacobianOrder = 3 * order - 1; break;
- case TYPE_HEX : jacobianOrder = 3 * order - 1; break;
+ 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;
default :
- Msg::Error("Unknown Jacobian function space for element type %d: reverting to TET_4", tag);
+ Msg::Error("Unknown Jacobian function space for element type %d", tag);
jacobianOrder = 0;
break;
}
- J.bezier = bezierBasis::find(polynomialBasis::getTag(parentType, jacobianOrder, false));
- polynomialBasis *F = (polynomialBasis*)BasisFactory::create(tag);
- generateGradShapes(J, J.bezier->points, F->monomials, F->coefficients);
+ jacType = polynomialBasis::getTag(parentType, jacobianOrder, false);
+ primJacType = polynomialBasis::getTag(parentType, primJacobianOrder, false);
}
+ J.bezier = bezierBasis::find(jacType);
+ const nodalBasis *mapBasis = BasisFactory::create(tag);
+ generateGradShapes(J, J.bezier->points, mapBasis);
+ const nodalBasis *primJacBasis = BasisFactory::create(primJacType);
+ generatePrimJac2Jac(J, primJacBasis);
return &J;
}
+
+
+
+void JacobianBasis::getSignedJacobian(MElement *el, fullVector<double> &jacobian) const
+{
+
+ const int numVertices = gradShapeMatX.size2(), numJacNodes = gradShapeMatX.size1();
+
+ switch (el->getDim()) {
+
+ case 1 : {
+ fullVector<double> nodesX(numVertices);
+ for (int i = 0; i < numVertices; i++) {
+ nodesX(i) = el->getShapeFunctionNode(i)->x();
+ }
+ 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();
+ }
+ fullMatrix<double> dxydX(numJacNodes,2), dxydY(numJacNodes,2);
+ gradShapeMatX.mult(nodesXY, dxydX);
+ gradShapeMatY.mult(nodesXY, 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();
+ }
+ fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3);
+ gradShapeMatX.mult(nodesXYZ, dxyzdX);
+ gradShapeMatY.mult(nodesXYZ, dxyzdY);
+ gradShapeMatZ.mult(nodesXYZ, dxyzdZ);
+ for (int i = 0; i < numJacNodes; i++) {
+ const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
+ const double &dxdY = dxyzdY(i,0), &dydY = dxyzdY(i,1), &dzdY = dxyzdY(i,2);
+ const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
+ jacobian(i) = dxdX*dydY*dzdZ + dxdY*dydZ*dzdX + dydX*dzdY*dxdZ
+ - dxdZ*dydY*dzdX - dxdY*dydX*dzdZ - dydZ*dzdY*dxdX;
+ }
+ break;
+ }
+
+ }
+
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index af8326a968153c54ea6a901bd7db55b7b4995da2..c7cf23c4785deb236ece3673edfdffac2937f93a 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -9,6 +9,9 @@
#include <map>
#include <vector>
#include "fullMatrix.h"
+
+class MElement;
+
class bezierBasis {
private :
static std::map<int, bezierBasis> _bbs;
@@ -29,10 +32,10 @@ class JacobianBasis {
static std::map<int, JacobianBasis> _fs;
public :
const bezierBasis *bezier;
- fullMatrix<double> gradShapeMatX;
- fullMatrix<double> gradShapeMatY;
- fullMatrix<double> gradShapeMatZ;
+ fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
+ fullMatrix<double> primJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
static const JacobianBasis *find(int);
+ void getSignedJacobian(MElement *el, fullVector<double> &jacobian) const;
};
#endif
diff --git a/Numeric/polynomialBasis.cpp b/Numeric/polynomialBasis.cpp
index 9334222f53aa690a80db01cf9e42737905645849..b50a8e78938e9cbcddc4213737fa047867da7c60 100644
--- a/Numeric/polynomialBasis.cpp
+++ b/Numeric/polynomialBasis.cpp
@@ -410,9 +410,10 @@ void polynomialBasis::f(const fullMatrix<double> &coord, fullMatrix<double> &sf)
sf.resize (coord.size1(), coefficients.size1());
for (int iPoint=0; iPoint< coord.size1(); iPoint++) {
evaluateMonomials(coord(iPoint,0), coord(iPoint,1), coord(iPoint,2), p);
- for (int i = 0; i < coefficients.size1(); i++)
- for (int j = 0; j < coefficients.size2(); j++)
- sf(iPoint,i) += coefficients(i, j) * p[j];
+ for (int i = 0; i < coefficients.size1(); i++) {
+ sf(iPoint,i) = 0.;
+ for (int j = 0; j < coefficients.size2(); j++) sf(iPoint,i) += coefficients(i, j) * p[j];
+ }
}
}
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index d577c875011df05742a858089cf3365cd1ea3efd..87036b090cb269966054012fee95fef3cdb3ded4 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -75,102 +75,6 @@ std::string GMSH_AnalyseCurvedMeshPlugin::getHelp() const
"Tolerance = R+ , << 1 : tolerance (relatively to J_min and J_max) used during the computation of J_min and J_max";
}
-// Miscellaneous method
-static void setJacobian(MElement *el, const JacobianBasis *jfs, fullVector<double> &jacobian)
-{
- int numVertices = el->getNumVertices();
- fullVector<double> nodesX(numVertices);
- fullVector<double> nodesY;
- fullVector<double> nodesZ;
- fullVector<double> interm1;
- fullVector<double> interm2;
-
- switch (el->getDim()) {
-
- case 1 :
- for (int i = 0; i < numVertices; i++) {
- nodesX(i) = el->getVertex(i)->x();
- }
- jfs->gradShapeMatX.mult(nodesX, jacobian);
- break;
-
- case 2 :
- nodesY.resize(numVertices);
- interm1.resize(jacobian.size());
- interm2.resize(jacobian.size());
-
- for (int i = 0; i < numVertices; i++) {
- nodesX(i) = el->getVertex(i)->x();
- nodesY(i) = el->getVertex(i)->y();
- }
-
- jfs->gradShapeMatX.mult(nodesX, jacobian);
- jfs->gradShapeMatY.mult(nodesY, interm2);
- jacobian.multTByT(interm2);
-
- jfs->gradShapeMatY.mult(nodesX, interm1);
- jfs->gradShapeMatX.mult(nodesY, interm2);
- interm1.multTByT(interm2);
-
- jacobian.axpy(interm1, -1);
- break;
-
- case 3 :
- nodesY.resize(numVertices);
- nodesZ.resize(numVertices);
- interm1.resize(jacobian.size());
- interm2.resize(jacobian.size());
-
- for (int i = 0; i < numVertices; i++) {
- nodesX(i) = el->getVertex(i)->x();
- nodesY(i) = el->getVertex(i)->y();
- nodesZ(i) = el->getVertex(i)->z();
- }
-
- jfs->gradShapeMatX.mult(nodesX, jacobian);
- jfs->gradShapeMatY.mult(nodesY, interm2);
- jacobian.multTByT(interm2);
- jfs->gradShapeMatZ.mult(nodesZ, interm2);
- jacobian.multTByT(interm2);
-
- jfs->gradShapeMatX.mult(nodesY, interm1);
- jfs->gradShapeMatY.mult(nodesZ, interm2);
- interm1.multTByT(interm2);
- jfs->gradShapeMatZ.mult(nodesX, interm2);
- interm1.multTByT(interm2);
- jacobian.axpy(interm1, 1);
-
- jfs->gradShapeMatX.mult(nodesZ, interm1);
- jfs->gradShapeMatY.mult(nodesX, interm2);
- interm1.multTByT(interm2);
- jfs->gradShapeMatZ.mult(nodesY, interm2);
- interm1.multTByT(interm2);
- jacobian.axpy(interm1, 1);
-
-
- jfs->gradShapeMatX.mult(nodesY, interm1);
- jfs->gradShapeMatY.mult(nodesX, interm2);
- interm1.multTByT(interm2);
- jfs->gradShapeMatZ.mult(nodesZ, interm2);
- interm1.multTByT(interm2);
- jacobian.axpy(interm1, -1);
-
- jfs->gradShapeMatX.mult(nodesZ, interm1);
- jfs->gradShapeMatY.mult(nodesY, interm2);
- interm1.multTByT(interm2);
- jfs->gradShapeMatZ.mult(nodesX, interm2);
- interm1.multTByT(interm2);
- jacobian.axpy(interm1, -1);
-
- jfs->gradShapeMatX.mult(nodesX, interm1);
- jfs->gradShapeMatY.mult(nodesZ, interm2);
- interm1.multTByT(interm2);
- jfs->gradShapeMatZ.mult(nodesY, interm2);
- interm1.multTByT(interm2);
- jacobian.axpy(interm1, -1);
- }
-}
-
static double sum(fullVector<double> &v)
{
double sum = .0;
@@ -334,8 +238,8 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
fullVector<double> jac1B(jfs1->bezier->points.size1(), numEl);
fullVector<double> jacBez, jacobian, jac1;
- setJacobian(el, jfs, jacobianB);
- setJacobian(el, jfs1, jac1B);
+ jfs->getSignedJacobian(el, jacobianB);
+ jfs1->getSignedJacobian(el, jac1B);
bb->matrixLag2Bez.mult(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
@@ -350,8 +254,8 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
jacobian.setAsProxy(jacobianB, k);
jac1.setAsProxy(jac1B, k);
#else
- setJacobian(el[k], jfs, jacobian);
- setJacobian(el[k], jfs1, jac1);
+ jfs->getSignedJacobian(el[k], jacobian);
+ jfs1->getSignedJacobian(el[k], jac1);
#endif
// AmJ : avgJ is not the average Jac for quad, prism or hex
@@ -539,8 +443,8 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
fullVector<double> jac1B(jfs1->bezier->points.size1(), numEl);
fullVector<double> jacBez, jacobian, jac1;
- setJacobian(el, jfs, jacobianB);
- setJacobian(el, jfs1, jac1B);
+ jfs->getSignedJacobian(el, jacobianB);
+ jfs1->getSignedJacobian(el, jac1B);
bb->matrixLag2Bez.mult(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
@@ -564,8 +468,8 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
jacobian.setAsProxy(jacobianB, k);
jac1.setAsProxy(jac1B, k);
#else
- setJacobian(el[k], jfs, jacobian);
- setJacobian(el[k], jfs1, jac1);
+ jfs->getSignedJacobian(el[k], jacobian);
+ jfs1->getSignedJacobian(el[k], jac1);
bb->matrixLag2Bez.mult(jacobian, jacBez);
#endif