diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index d2ea64588dfc226611ac108f1a6856915bb3a290..5e28c6fb00bb17399966e69251524d430dffd467 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -91,12 +91,12 @@ const MetricBasis* BasisFactory::getMetricBasis(int tag)
return M;
}
-const CondNumBasis* BasisFactory::getCondNumBasis(int tag)
+const CondNumBasis* BasisFactory::getCondNumBasis(int tag, int cnOrder)
{
std::map<int, CondNumBasis*>::const_iterator it = cs.find(tag);
if (it != cs.end()) return it->second;
- CondNumBasis* M = new CondNumBasis(tag);
+ CondNumBasis* M = new CondNumBasis(tag, cnOrder);
cs.insert(std::make_pair(tag, M));
return M;
}
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index ad408f48b7397b3819fc02078dca533e7188c4cf..d5c509c1640cd9e2dfd19738d84228d39acb6ad2 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -54,7 +54,7 @@ class BasisFactory
static const MetricBasis* getMetricBasis(int tag);
// Condition number
- static const CondNumBasis* getCondNumBasis(int tag);
+ static const CondNumBasis* getCondNumBasis(int tag, int cnOrder = -1);
// Gradients
static const GradientBasis* getGradientBasis(FuncSpaceData);
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index d67dcad36f4b4a2e2decb8fd059982fa1d46ebf5..c443f09fea2977f43992be7878210585232c8f89 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -12,6 +12,51 @@
namespace {
+template<class T>
+void calcMapFromIdealElement(int type, T &gSMatX, T &gSMatY, T &gSMatZ)
+{
+ // 2D scaling
+ switch(type) {
+ case TYPE_QUA: // Quad, hex, pyramid -> square with side of length 1
+ case TYPE_HEX:
+ case TYPE_PYR: {
+ gSMatX.scale(2.);
+ gSMatY.scale(2.);
+ break;
+ }
+ default: { // Tri, tet, prism: equilateral tri with side of length 1
+ static const double cTri[2] = {-1./std::sqrt(3), 2./std::sqrt(3)};
+ gSMatY.scale(cTri[1]);
+ gSMatY.axpy(gSMatX, cTri[0]);
+ break;
+ }
+ }
+
+ // 3D scaling
+ switch(type) {
+ case TYPE_HEX: // Hex, prism -> side of length 1 in z
+ case TYPE_PRI: {
+ gSMatZ.scale(2.);
+ break;
+ }
+ case TYPE_PYR: { // Pyramid -> height sqrt(2)/2
+ static const double cPyr = 1./sqrt(2.);
+ gSMatZ.scale(cPyr);
+ break;
+ }
+ case TYPE_TET: // Tet: take into account (x, y) scaling to obtain regular tet
+ {
+ static const double cTet[3] = {-3./2/std::sqrt(6),
+ -1./2/std::sqrt(2),
+ std::sqrt(1.5)};
+ gSMatZ.scale(cTet[2]);
+ gSMatZ.axpy(gSMatX, cTet[0]);
+ gSMatZ.axpy(gSMatY, cTet[1]);
+ break;
+ }
+ }
+}
+
// Compute the determinant of a 3x3 matrix
inline double calcDet3D(double M11, double M12, double M13,
double M21, double M22, double M23,
@@ -134,8 +179,8 @@ GradientBasis::GradientBasis(FuncSpaceData data)
gradShapeIdealMatX = gradShapeMatX;
gradShapeIdealMatY = gradShapeMatY;
gradShapeIdealMatZ = gradShapeMatZ;
- mapFromIdealElement(_data.elementType(), &gradShapeIdealMatX,
- &gradShapeIdealMatY, &gradShapeIdealMatZ);
+ mapFromIdealElement(_data.elementType(), gradShapeIdealMatX,
+ gradShapeIdealMatY, gradShapeIdealMatZ);
}
void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
@@ -159,56 +204,25 @@ void GradientBasis::getIdealGradientsFromNodes(const fullMatrix<double> &nodes,
}
void GradientBasis::mapFromIdealElement(int type,
- fullMatrix<double> *dxyzdX,
- fullMatrix<double> *dxyzdY,
- fullMatrix<double> *dxyzdZ)
+ fullMatrix<double> &gSMatX,
+ fullMatrix<double> &gSMatY,
+ fullMatrix<double> &gSMatZ)
{
- // 2D scaling
- switch(type) {
- case TYPE_QUA: // Quad, hex, pyramid -> square with side of length 1
- case TYPE_HEX:
- case TYPE_PYR: {
- dxyzdX->scale(2.);
- dxyzdY->scale(2.);
- break;
- }
- default: { // Tri, tet, prism: equilateral tri with side of length 1
- static const double cTri[2] = {-1./std::sqrt(3), 2./std::sqrt(3)};
- dxyzdY->scale(cTri[1]);
- dxyzdY->axpy(*dxyzdX, cTri[0]);
- break;
- }
- }
+ calcMapFromIdealElement(type, gSMatX, gSMatY, gSMatZ);
+}
- // 3D scaling
- switch(type) {
- case TYPE_HEX: // Hex, prism -> side of length 1 in z
- case TYPE_PRI: {
- dxyzdZ->scale(2.);
- break;
- }
- case TYPE_PYR: { // Pyramid -> height sqrt(2)/2
- static const double cPyr = 1./sqrt(2.);
- dxyzdZ->scale(cPyr);
- break;
- }
- case TYPE_TET: // Tet: take into account (x, y) scaling to obtain regular tet
- {
- static const double cTet[3] = {-3./2/std::sqrt(6),
- -1./2/std::sqrt(2),
- std::sqrt(1.5)};
- dxyzdZ->scale(cTet[2]);
- dxyzdZ->axpy(*dxyzdX, cTet[0]);
- dxyzdZ->axpy(*dxyzdY, cTet[1]);
- break;
- }
- }
+void GradientBasis::mapFromIdealElement(int type,
+ fullVector<double> &gSVecX,
+ fullVector<double> &gSVecY,
+ fullVector<double> &gSVecZ)
+{
+ calcMapFromIdealElement(type, gSVecX, gSVecY, gSVecZ);
}
void GradientBasis::mapFromIdealElement(int type, double jac[3][3])
{
fullMatrix<double> dxyzdX(jac[0], 1, 3), dxyzdY(jac[1], 1, 3), dxyzdZ(jac[2], 1, 3);
- mapFromIdealElement(type, &dxyzdX, &dxyzdY, &dxyzdZ);
+ mapFromIdealElement(type, dxyzdX, dxyzdY, dxyzdZ);
}
JacobianBasis::JacobianBasis(FuncSpaceData data)
@@ -252,15 +266,21 @@ JacobianBasis::JacobianBasis(FuncSpaceData data)
double (*barDPsi)[3] = new double[numPrimMapNodes][3];
primMapBasis->df(xBar, yBar, zBar, barDPsi);
- primGradShapeBarycenterX.resize(numPrimMapNodes);
- primGradShapeBarycenterY.resize(numPrimMapNodes);
- primGradShapeBarycenterZ.resize(numPrimMapNodes);
+ primGradShapeBaryX.resize(numPrimMapNodes);
+ primGradShapeBaryY.resize(numPrimMapNodes);
+ primGradShapeBaryZ.resize(numPrimMapNodes);
for (int j = 0; j < numPrimMapNodes; j++) {
- primGradShapeBarycenterX(j) = barDPsi[j][0];
- primGradShapeBarycenterY(j) = barDPsi[j][1];
- primGradShapeBarycenterZ(j) = barDPsi[j][2];
+ primGradShapeBaryX(j) = barDPsi[j][0];
+ primGradShapeBaryY(j) = barDPsi[j][1];
+ primGradShapeBaryZ(j) = barDPsi[j][2];
}
+ primIdealGradShapeBaryX = primGradShapeBaryX;
+ primIdealGradShapeBaryY = primGradShapeBaryY;
+ primIdealGradShapeBaryZ = primGradShapeBaryZ;
+ _gradBasis->mapFromIdealElement(primIdealGradShapeBaryX, primIdealGradShapeBaryY,
+ primIdealGradShapeBaryZ);
+
delete[] barDPsi;
// Compute "fast" Jacobian evaluation matrices (at 1st order nodes + barycenter)
@@ -298,14 +318,15 @@ const bezierBasis* JacobianBasis::getBezier() const
}
// Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value)
-double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
+double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ,
+ fullMatrix<double> &result) const
{
fullVector<double> dxyzdXbar(3);
for (int j = 0; j < numPrimMapNodes; j++) {
- dxyzdXbar(0) += primGradShapeBarycenterX(j) * nodesXYZ(j, 0);
- dxyzdXbar(1) += primGradShapeBarycenterX(j) * nodesXYZ(j, 1);
- dxyzdXbar(2) += primGradShapeBarycenterX(j) * nodesXYZ(j, 2);
+ dxyzdXbar(0) += primGradShapeBaryX(j) * nodesXYZ(j, 0);
+ dxyzdXbar(1) += primGradShapeBaryX(j) * nodesXYZ(j, 1);
+ dxyzdXbar(2) += primGradShapeBaryX(j) * nodesXYZ(j, 2);
}
if((fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(1)) && fabs(dxyzdXbar(0)) >= fabs(dxyzdXbar(2))) ||
@@ -329,19 +350,21 @@ double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullM
}
// Computes (unit) normal to straight surface element at barycenter (with norm as return value)
-double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result, bool ideal) const
+double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ,
+ fullMatrix<double> &result, bool ideal) const
{
+ const fullVector<double> &gSX = ideal ? primIdealGradShapeBaryX : primGradShapeBaryX,
+ &gSY = ideal ? primIdealGradShapeBaryY : primGradShapeBaryY;
fullMatrix<double> dxyzdX(1, 3), dxyzdY(1, 3);
for (int j=0; j<numPrimMapNodes; j++) {
- dxyzdX(0, 0) += primGradShapeBarycenterX(j)*nodesXYZ(j,0);
- dxyzdX(0, 1) += primGradShapeBarycenterX(j)*nodesXYZ(j,1);
- dxyzdX(0, 2) += primGradShapeBarycenterX(j)*nodesXYZ(j,2);
- dxyzdY(0, 0) += primGradShapeBarycenterY(j)*nodesXYZ(j,0);
- dxyzdY(0, 1) += primGradShapeBarycenterY(j)*nodesXYZ(j,1);
- dxyzdY(0, 2) += primGradShapeBarycenterY(j)*nodesXYZ(j,2);
+ dxyzdX(0, 0) += gSX(j)*nodesXYZ(j,0);
+ dxyzdX(0, 1) += gSX(j)*nodesXYZ(j,1);
+ dxyzdX(0, 2) += gSX(j)*nodesXYZ(j,2);
+ dxyzdY(0, 0) += gSY(j)*nodesXYZ(j,0);
+ dxyzdY(0, 1) += gSY(j)*nodesXYZ(j,1);
+ dxyzdY(0, 2) += gSY(j)*nodesXYZ(j,2);
}
- if (ideal) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, 0);
result(0, 2) = dxyzdX(0, 0) * dxyzdY(0, 1) - dxyzdX(0, 1) * dxyzdY(0, 0);
result(0, 1) = -dxyzdX(0, 0) * dxyzdY(0, 2) + dxyzdX(0, 2) * dxyzdY(0, 0);
result(0, 0) = dxyzdX(0, 1) * dxyzdY(0, 2) - dxyzdX(0, 2) * dxyzdY(0, 1);
@@ -355,20 +378,22 @@ double JacobianBasis::getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMa
// Returns absolute value of Jacobian of straight volume element at barycenter
double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ, bool ideal) const
{
+ const fullVector<double> &gSX = ideal ? primIdealGradShapeBaryX : primGradShapeBaryX,
+ &gSY = ideal ? primIdealGradShapeBaryY : primGradShapeBaryY,
+ &gSZ = ideal ? primIdealGradShapeBaryZ : primGradShapeBaryZ;
fullMatrix<double> dxyzdX(1, 3), dxyzdY(1, 3), dxyzdZ(1, 3);
for (int j=0; j<numPrimMapNodes; j++) {
- dxyzdX(0, 0) += primGradShapeBarycenterX(j)*nodesXYZ(j,0);
- dxyzdX(0, 1) += primGradShapeBarycenterX(j)*nodesXYZ(j,1);
- dxyzdX(0, 2) += primGradShapeBarycenterX(j)*nodesXYZ(j,2);
- dxyzdY(0, 0) += primGradShapeBarycenterY(j)*nodesXYZ(j,0);
- dxyzdY(0, 1) += primGradShapeBarycenterY(j)*nodesXYZ(j,1);
- dxyzdY(0, 2) += primGradShapeBarycenterY(j)*nodesXYZ(j,2);
- dxyzdZ(0, 0) += primGradShapeBarycenterZ(j)*nodesXYZ(j,0);
- dxyzdZ(0, 1) += primGradShapeBarycenterZ(j)*nodesXYZ(j,1);
- dxyzdZ(0, 2) += primGradShapeBarycenterZ(j)*nodesXYZ(j,2);
+ dxyzdX(0, 0) += gSX(j)*nodesXYZ(j,0);
+ dxyzdX(0, 1) += gSX(j)*nodesXYZ(j,1);
+ dxyzdX(0, 2) += gSX(j)*nodesXYZ(j,2);
+ dxyzdY(0, 0) += gSY(j)*nodesXYZ(j,0);
+ dxyzdY(0, 1) += gSY(j)*nodesXYZ(j,1);
+ dxyzdY(0, 2) += gSY(j)*nodesXYZ(j,2);
+ dxyzdZ(0, 0) += gSZ(j)*nodesXYZ(j,0);
+ dxyzdZ(0, 1) += gSZ(j)*nodesXYZ(j,1);
+ dxyzdZ(0, 2) += gSZ(j)*nodesXYZ(j,2);
}
- if (ideal) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, &dxyzdZ);
double dJ = fabs(calcDet3D(dxyzdX(0, 0), dxyzdY(0, 0), dxyzdZ(0, 0),
dxyzdX(0, 1), dxyzdY(0, 1), dxyzdZ(0, 1),
dxyzdX(0, 2), dxyzdY(0, 2), dxyzdZ(0, 2)));
@@ -377,10 +402,14 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ, bool idea
// Calculate (signed, possibly scaled) Jacobian for one element, with normal vectors to straight element
// for regularization. Evaluation points depend on the given matrices for shape function gradients.
-template<bool scaling, bool ideal>
-inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
+template<bool scaling>
+inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ bool idealNorm,
+ fullVector<double> &jacobian) const
{
switch (_dim) {
@@ -411,7 +440,7 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
case 2 : {
fullMatrix<double> normal(1,3);
- const double invScale = getPrimNormal2D(nodesXYZ,normal);
+ const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
if (scaling) {
const double scale = 1./invScale;
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
@@ -419,7 +448,6 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
- if (ideal) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, NULL);
for (int i = 0; i < nJacNodes; 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);
@@ -437,7 +465,6 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
gSMatZ.mult(nodesXYZ, dxyzdZ);
- if (ideal) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, &dxyzdZ);
for (int i = 0; i < nJacNodes; 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);
@@ -460,39 +487,43 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
// Calculate signed Jacobian for one element, with normal vectors to straight element for
// regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
-{
- getJacobianGeneral<false, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
-}
-
-// Calculate signed Jacobian (from ideal) for one element, with normal vectors to straight element for
-// regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getSignedIdealJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
+void JacobianBasis::getSignedJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ bool idealNorm,
+ fullVector<double> &jacobian) const
{
- getJacobianGeneral<false, true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
+ getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, idealNorm, jacobian);
}
// Calculate (signed) scaled Jacobian for one element, with normal vectors to straight element
// for regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
+void JacobianBasis::getScaledJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ bool idealNorm,
+ fullVector<double> &jacobian) const
{
- getJacobianGeneral<true, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
+ getJacobianGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, idealNorm, jacobian);
}
// Calculate (signed, possibly scaled) Jacobian for several elements, with normal vectors to straight
// elements for regularization. Evaluation points depend on the given matrices for shape function gradients.
// TODO: Optimize and test 1D & 2D
-template<bool scaling, bool ideal>
-inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
+template<bool scaling>
+inline void JacobianBasis::getJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ bool idealNorm,
+ fullMatrix<double> &jacobian) const
{
switch (_dim) {
@@ -536,11 +567,6 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
fullMatrix<double> dxdY(nJacNodes,numEl), dydY(nJacNodes,numEl), dzdY(nJacNodes,numEl);
gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
- if (ideal) {
- _gradBasis->mapFromIdealElement(&dxdX, &dxdY, NULL);
- _gradBasis->mapFromIdealElement(&dydX, &dydY, NULL);
- _gradBasis->mapFromIdealElement(&dzdX, &dzdY, NULL);
- }
for (int iEl = 0; iEl < numEl; iEl++) {
fullMatrix<double> nodesXYZ(numPrimMapNodes,3);
for (int i = 0; i < numPrimMapNodes; i++) {
@@ -549,7 +575,7 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normal(1,3);
- const double invScale = getPrimNormal2D(nodesXYZ,normal);
+ const double invScale = getPrimNormal2D(nodesXYZ, normal, idealNorm);
if (scaling) {
const double scale = 1./invScale;
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
@@ -571,11 +597,6 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
gSMatZ.mult(nodesX, dxdZ); gSMatZ.mult(nodesY, dydZ); gSMatZ.mult(nodesZ, dzdZ);
- if (ideal) {
- _gradBasis->mapFromIdealElement(&dxdX, &dxdY, &dxdZ);
- _gradBasis->mapFromIdealElement(&dydX, &dydY, &dydZ);
- _gradBasis->mapFromIdealElement(&dzdX, &dzdY, &dzdZ);
- }
for (int iEl = 0; iEl < numEl; iEl++) {
for (int i = 0; i < nJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl), dxdY(i,iEl), dxdZ(i,iEl),
@@ -601,37 +622,38 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
// Calculate signed Jacobian for several elements, with normal vectors to straight elements for
// regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
-{
- getJacobianGeneral<false, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
-}
-
-// Calculate signed Jacobian (from ideal) for several elements, with normal vectors to straight elements
-// for regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getSignedIdealJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
+void JacobianBasis::getSignedJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ bool idealNorm,
+ fullMatrix<double> &jacobian) const
{
- getJacobianGeneral<false, true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
+ getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ,
+ nodesX, nodesY, nodesZ, idealNorm, jacobian);
}
// Calculate (signed) scaled Jacobian for several elements, with normal vectors to straight elements
// for regularization. Evaluation points depend on the given matrices for shape function gradients.
-void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
+void JacobianBasis::getScaledJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ bool idealNorm,
+ fullMatrix<double> &jacobian) const
{
- getJacobianGeneral<true, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
+ getJacobianGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ,
+ nodesX, nodesY, nodesZ, idealNorm, jacobian);
}
// Calculate (signed) Jacobian and its gradients for one element, with normal vectors to straight element
// for regularization. Evaluation points depend on the given matrices for shape function gradients.
-template<bool ideal>
inline void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY,
@@ -674,7 +696,6 @@ inline void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
- if (ideal) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, NULL);
for (int i = 0; i < nJacNodes; 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);
@@ -693,7 +714,6 @@ inline void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
gSMatZ.mult(nodesXYZ, dxyzdZ);
- if (ideal) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, &dxyzdZ);
for (int i = 0; i < nJacNodes; 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);
@@ -711,17 +731,6 @@ inline void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
}
-void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
- const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY,
- const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals,
- fullMatrix<double> &JDJ) const
-{
- getSignedJacAndGradientsGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, normals, JDJ);
-}
-
void JacobianBasis::getSignedIdealJacAndGradientsGeneral(int nJacNodes,
const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY,
@@ -730,7 +739,7 @@ void JacobianBasis::getSignedIdealJacAndGradientsGeneral(int nJacNodes,
const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const
{
- getSignedJacAndGradientsGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, normals, JDJ);
+ getSignedJacAndGradientsGeneral(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, normals, JDJ);
}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 978e0cd1e43dd739f080e23f183441d295d95ac1..7bfbf853757ea510a186319f888c3b3be66433ec 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -32,15 +32,24 @@ public:
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
fullMatrix<double> *dxyzdZ) const;
- void mapFromIdealElement(fullMatrix<double> *dxyzdX,
- fullMatrix<double> *dxyzdY,
- fullMatrix<double> *dxyzdZ) const {
+ void mapFromIdealElement(fullMatrix<double> &dxyzdX,
+ fullMatrix<double> &dxyzdY,
+ fullMatrix<double> &dxyzdZ) const {
+ GradientBasis::mapFromIdealElement(_data.elementType(), dxyzdX, dxyzdY, dxyzdZ);
+ }
+ void mapFromIdealElement(fullVector<double> &dxyzdX,
+ fullVector<double> &dxyzdY,
+ fullVector<double> &dxyzdZ) const {
GradientBasis::mapFromIdealElement(_data.elementType(), dxyzdX, dxyzdY, dxyzdZ);
}
static void mapFromIdealElement(int type,
- fullMatrix<double> *dxyzdX,
- fullMatrix<double> *dxyzdY,
- fullMatrix<double> *dxyzdZ);
+ fullMatrix<double> &gSMatX,
+ fullMatrix<double> &gSMatY,
+ fullMatrix<double> &gSMatZ);
+ static void mapFromIdealElement(int type,
+ fullVector<double> &gSVecX,
+ fullVector<double> &gSVecY,
+ fullVector<double> &gSVecZ);
static void mapFromIdealElement(int type, double jac[3][3]);
};
@@ -53,8 +62,10 @@ private:
const int _dim;
fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
- fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
- fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
+ fullVector<double> primGradShapeBaryX, primGradShapeBaryY, primGradShapeBaryZ;
+ fullVector<double> primIdealGradShapeBaryX, primIdealGradShapeBaryY,
+ primIdealGradShapeBaryZ;
+ fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
int numJacNodes, numPrimJacNodes;
int numMapNodes, numPrimMapNodes;
@@ -73,67 +84,105 @@ public:
const bezierBasis* getBezier() const;
// Jacobian evaluation methods
- double getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
- double getPrimNormal2D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result, bool ideal=false) const;
+ double getPrimNormals1D(const fullMatrix<double> &nodesXYZ,
+ fullMatrix<double> &result) const;
+ double getPrimNormal2D(const fullMatrix<double> &nodesXYZ,
+ fullMatrix<double> &result, bool ideal=false) const;
double getPrimJac3D(const fullMatrix<double> &nodesXYZ, bool ideal=false) const;
inline void getSignedJacAndGradients(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &JDJ) const {
getSignedJacAndGradientsGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ, nodesXYZ, normals, JDJ);
+ _gradBasis->gradShapeMatY,
+ _gradBasis->gradShapeMatZ,
+ nodesXYZ, normals, JDJ);
}
inline void getSignedJacAndGradientsFast(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
- getSignedJacAndGradientsGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
- gradShapeMatZFast,nodesXYZ,normals,JDJ);
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &JDJ) const {
+ getSignedJacAndGradientsGeneral(numJacNodesFast, gradShapeMatXFast,
+ gradShapeMatYFast, gradShapeMatZFast,
+ nodesXYZ, normals, JDJ);
}
inline void getSignedIdealJacAndGradients(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
- getSignedIdealJacAndGradientsGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ, nodesXYZ, normals, JDJ);
- }
- inline void getSignedIdealJacAndGradientsFast(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
- getSignedIdealJacAndGradientsGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
- gradShapeMatZFast,nodesXYZ,normals,JDJ);
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &JDJ) const {
+ getSignedJacAndGradientsGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
+ _gradBasis->gradShapeIdealMatY,
+ _gradBasis->gradShapeIdealMatZ,
+ nodesXYZ, normals, JDJ);
}
+// inline void getSignedIdealJacAndGradientsFast(const fullMatrix<double> &nodesXYZ,
+// const fullMatrix<double> &normals, fullMatrix<double> &JDJ) const {
+// getSignedIdealJacAndGradientsGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
+// gradShapeMatZFast,nodesXYZ,normals,JDJ);
+// }
void getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &nodesXYZStraight,
- fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const;
- inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
+ fullVector<double> &lambdaJ,
+ fullMatrix<double> &gradLambdaJ) const;
+ inline void getSignedJacobian(const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesXYZ,jacobian);
+ _gradBasis->gradShapeMatY,
+ _gradBasis->gradShapeMatZ,
+ nodesXYZ, false, jacobian);
}
- inline void getSignedJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
+ inline void getSignedJacobian(const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ fullMatrix<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
+ _gradBasis->gradShapeMatY,
+ _gradBasis->gradShapeMatZ,
+ nodesX, nodesY, nodesZ, false, jacobian);
}
- inline void getSignedIdealJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
- getSignedIdealJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesXYZ,jacobian);
+ inline void getSignedIdealJacobian(const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &jacobian) const {
+ getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
+ _gradBasis->gradShapeIdealMatY,
+ _gradBasis->gradShapeIdealMatZ,
+ nodesXYZ, true, jacobian);
}
- inline void getSignedIdealJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
- getSignedIdealJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
+ inline void getSignedIdealJacobian(const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ fullMatrix<double> &jacobian) const {
+ getSignedJacobianGeneral(numJacNodes, _gradBasis->gradShapeIdealMatX,
+ _gradBasis->gradShapeIdealMatY,
+ _gradBasis->gradShapeIdealMatZ,
+ nodesX, nodesY, nodesZ, true, jacobian);
}
- inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
+ inline void getScaledJacobian(const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesXYZ,jacobian);
+ _gradBasis->gradShapeMatY,
+ _gradBasis->gradShapeMatZ,
+ nodesXYZ, false, jacobian);
}
- inline void getScaledJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
+ inline void getScaledJacobian(const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ fullMatrix<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
- }
- inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
- getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
+ _gradBasis->gradShapeMatY,
+ _gradBasis->gradShapeMatZ,
+ nodesX, nodesY, nodesZ, false, jacobian);
}
- inline void getSignedIdealJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
- getSignedIdealJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
+ inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &jacobian) const {
+ getSignedJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
+ gradShapeMatYFast, gradShapeMatZFast,
+ nodesXYZ, false, jacobian);
}
- inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
- getScaledJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
+// inline void getSignedIdealJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
+// getSignedIdealJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
+// }
+ inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &jacobian) const {
+ getScaledJacobianGeneral(numJacNodesFast, gradShapeMatXFast,
+ gradShapeMatYFast, gradShapeMatZFast,
+ nodesXYZ, false, jacobian);
}
//
inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
@@ -156,49 +205,64 @@ public:
private :
- template<bool scaling, bool ideal>
- void getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
- template<bool scaling, bool ideal>
- void getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
- void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
- void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
- void getSignedIdealJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
- void getSignedIdealJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
- void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
- void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
- const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
-
- template<bool ideal>
- void getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
- fullMatrix<double> &JDJ) const;
- void getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
+ template<bool scaling>
+ void getJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ bool idealNorm, fullVector<double> &jacobian) const;
+ template<bool scaling>
+ void getJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ bool idealNorm, fullMatrix<double> &jacobian) const;
+ void getSignedJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ bool idealNorm, fullVector<double> &jacobian) const;
+ void getSignedJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ bool idealNorm, fullMatrix<double> &jacobian) const;
+ void getScaledJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ bool idealNorm, fullVector<double> &jacobian) const;
+ void getScaledJacobianGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesX,
+ const fullMatrix<double> &nodesY,
+ const fullMatrix<double> &nodesZ,
+ bool idealNorm, fullMatrix<double> &jacobian) const;
+
+ void getSignedJacAndGradientsGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const;
- void getSignedIdealJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
+ void getSignedIdealJacAndGradientsGeneral(int nJacNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const;
};
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 4f7837502161bd536ea709f09e6946b285b22d19..eb4959eec41d8a1bfb7058f82b23c0b0567125b0 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -1121,8 +1121,10 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
case 2 :
{
fullMatrix<double> dxydX(nSampPnts,3), dxydY(nSampPnts,3);
- gradients->getGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
- if (ideal) gradients->mapFromIdealElement(&dxydX, &dxydY, NULL);
+ if (ideal)
+ gradients->getIdealGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
+ else
+ gradients->getGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
coeff.resize(nSampPnts, 3);
for (int i = 0; i < nSampPnts; i++) {
@@ -1149,8 +1151,10 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
case 3 :
{
fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
- gradients->getGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
- if (ideal) gradients->mapFromIdealElement(&dxyzdX, &dxyzdY, &dxyzdZ);
+ if (ideal)
+ gradients->getIdealGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
+ else
+ gradients->getGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
coeff.resize(nSampPnts, 7);
for (int i = 0; i < nSampPnts; i++) {