diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index a9fc13c6bd4b4db88c5d8c82a7108417ffc8bec5..dcc9c52a2ebb2c603f93284d02a6b1f03f01225f 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -22,7 +22,7 @@
#include "GEntity.h"
#include "StringUtils.h"
#include "Numeric.h"
-#include "MetricBasis.h"
+#include "CondNumBasis.h"
#include "Context.h"
#define SQU(a) ((a)*(a))
@@ -322,6 +322,52 @@ void MElement::idealJacRange(double &jmin, double &jmax, GEntity *ge)
#endif
}
+void MElement::invCondNumRange(double &iCNMin, double &iCNMax, GEntity *ge)
+{
+ iCNMin = iCNMax = 1.0;
+#if defined(HAVE_MESH)
+ const CondNumBasis *cnb = BasisFactory::getCondNumBasis(getTypeForMSH());
+ const int numCNNodes = cnb->getNumCondNumNodes();
+ fullMatrix<double> nodesXYZ(cnb->getNumMapNodes(), 3), normals;
+ getNodesCoord(nodesXYZ);
+ if (getDim() == 2.) {
+ SVector3 nVec = getFace(0).normal();
+ normals.resize(1, 3);
+ normals(0, 0) = nVec[0]; normals(0, 1) = nVec[1]; normals(0, 2) = nVec[2];
+ }
+ if (ge && (ge->dim() == 2) && ge->haveParametrization()) {
+ // If parametrized surface entity provided...
+ SVector3 geoNorm(0., 0., 0.);
+ // ... correct Jacobian sign with geometrical normal
+ for (int i=0; i<getNumPrimaryVertices(); i++) {
+ const MVertex *vert = getVertex(i);
+ if (vert->onWhat() == ge) {
+ double u, v;
+ vert->getParameter(0, u);
+ vert->getParameter(1, v);
+ geoNorm += ((GFace*)ge)->normal(SPoint2(u, v));
+ }
+ }
+ if (geoNorm.normSq() == 0.) {
+ // If no vertex on surface or average is zero, take normal at barycenter
+ SPoint2 param = ((GFace*)ge)->parFromPoint(barycenter(true), false);
+ geoNorm = ((GFace*)ge)->normal(param);
+ }
+ const double dp = geoNorm(0) * normals(0, 0) + geoNorm(1) * normals(0, 1)
+ + geoNorm(2) * normals(0, 2);
+ if (dp < 0.) {
+ normals(0, 0) = -normals(0, 0);
+ normals(0, 1) = -normals(0, 1);
+ normals(0, 2) = -normals(0, 2);
+ }
+ }
+ fullVector<double> invCondNum(numCNNodes);
+ cnb->getSignedInvCondNum(nodesXYZ, normals, invCondNum);
+ iCNMin = *std::min_element(invCondNum.getDataPtr(), invCondNum.getDataPtr()+numCNNodes);
+ iCNMax = *std::max_element(invCondNum.getDataPtr(), invCondNum.getDataPtr()+numCNNodes);
+#endif
+}
+
void MElement::getNode(int num, double &u, double &v, double &w) const
{
// only for MElements that don't have a lookup table for this
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 5b92db8693879d5d932fa0b1347b3211e6a04636..1ce7cfdeb6e6a5dfa6bc8bbfa45234f14f3b2da2 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -17,6 +17,7 @@
#include "nodalBasis.h"
#include "polynomialBasis.h"
#include "JacobianBasis.h"
+#include "MetricBasis.h"
#include "GaussIntegration.h"
class GModel;
@@ -204,6 +205,7 @@ class MElement
virtual double angleShapeMeasure() { return 1.0; }
virtual void scaledJacRange(double &jmin, double &jmax, GEntity *ge = 0) const;
virtual void idealJacRange(double &jmin, double &jmax, GEntity *ge = 0);
+ virtual void invCondNumRange(double &iCNMin, double &iCNMax, GEntity *ge = 0);
virtual double metricShapeMeasure();
virtual double metricShapeMeasure2();
diff --git a/Numeric/BasisFactory.cpp b/Numeric/BasisFactory.cpp
index d5678f959dd998ca1ce1148913954ee97269fc0e..d2ea64588dfc226611ac108f1a6856915bb3a290 100644
--- a/Numeric/BasisFactory.cpp
+++ b/Numeric/BasisFactory.cpp
@@ -9,11 +9,13 @@
#include "pyramidalBasis.h"
#include "miniBasis.h"
#include "MetricBasis.h"
+#include "CondNumBasis.h"
#include <map>
#include <cstddef>
std::map<int, nodalBasis*> BasisFactory::fs;
std::map<int, MetricBasis*> BasisFactory::ms;
+std::map<int, CondNumBasis*> BasisFactory::cs;
std::map<FuncSpaceData, JacobianBasis*> BasisFactory::js;
std::map<FuncSpaceData, bezierBasis*> BasisFactory::bs;
std::map<FuncSpaceData, GradientBasis*> BasisFactory::gs;
@@ -89,6 +91,16 @@ const MetricBasis* BasisFactory::getMetricBasis(int tag)
return M;
}
+const CondNumBasis* BasisFactory::getCondNumBasis(int tag)
+{
+ std::map<int, CondNumBasis*>::const_iterator it = cs.find(tag);
+ if (it != cs.end()) return it->second;
+
+ CondNumBasis* M = new CondNumBasis(tag);
+ cs.insert(std::make_pair(tag, M));
+ return M;
+}
+
const GradientBasis* BasisFactory::getGradientBasis(FuncSpaceData data)
{
std::map<FuncSpaceData, GradientBasis*>::const_iterator it = gs.find(data);
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 9d5ac2accac357c525db645a07f8a4c59a24d80e..ad408f48b7397b3819fc02078dca533e7188c4cf 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -12,12 +12,14 @@ class nodalBasis;
class MetricBasis;
class GradientBasis;
class bezierBasis;
+class CondNumBasis;
class BasisFactory
{
private:
static std::map<int, nodalBasis*> fs;
static std::map<int, MetricBasis*> ms;
+ static std::map<int, CondNumBasis*> cs;
static std::map<FuncSpaceData, JacobianBasis*> js;
static std::map<FuncSpaceData, bezierBasis*> bs;
static std::map<FuncSpaceData, GradientBasis*> gs;
@@ -51,6 +53,9 @@ class BasisFactory
// Metric
static const MetricBasis* getMetricBasis(int tag);
+ // Condition number
+ static const CondNumBasis* getCondNumBasis(int tag);
+
// Gradients
static const GradientBasis* getGradientBasis(FuncSpaceData);
static const GradientBasis* getGradientBasis(int tag, int order) {
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index a865d9416e241c08fc6e9a934147d0bf846ba9dd..0cc2843650854b7a223a489eb8beed6e2bd09e43 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -19,6 +19,7 @@ set(SRC
bezierBasis.cpp
JacobianBasis.cpp
MetricBasis.cpp
+ CondNumBasis.cpp
pointsGenerators.cpp
ElementType.cpp
GaussIntegration.cpp
diff --git a/Numeric/CondNumBasis.cpp b/Numeric/CondNumBasis.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ce48ccef9902f8062aedd6edcca2274cf337fa9d
--- /dev/null
+++ b/Numeric/CondNumBasis.cpp
@@ -0,0 +1,584 @@
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#include "CondNumBasis.h"
+
+#include "GmshDefines.h"
+#include "GmshMessage.h"
+
+#include <vector>
+#include "polynomialBasis.h"
+#include "pyramidalBasis.h"
+#include "pointsGenerators.h"
+#include "BasisFactory.h"
+#include "Numeric.h"
+
+
+namespace {
+
+
+// Compute the determinant of a 3x3 matrix
+inline double calcDet3x3(double M11, double M12, double M13,
+ double M21, double M22, double M23,
+ double M31, double M32, double M33)
+{
+ return M11 * (M22*M33 - M23*M32)
+ - M12 * (M21*M33 - M23*M31)
+ + M13 * (M21*M32 - M22*M31);
+}
+
+
+// Compute the squared Frobenius norm of the inverse of a matrix
+template<bool sign>
+inline double calcInvCondNum2D(double dxdX, double dxdY,
+ double dydX, double dydY,
+ double dzdX, double dzdY,
+ double nx, double ny, double nz)
+{
+ static const double Eps = 1e-10;
+ const double dxdXSq = dxdX*dxdX, dydXSq = dydX*dydX, dzdXSq = dzdX*dzdX;
+ const double dxdYSq = dxdY*dxdY, dydYSq = dydY*dydY, dzdYSq = dzdY*dzdY;
+ const double Dx = dxdXSq - dxdYSq, Dy = dydXSq - dydYSq;
+ const double Cx = dxdX * dxdY, Cy = dydX * dydY;
+ const double S1 = dzdYSq * dzdYSq;
+ const double S2 = (dzdXSq - Dy - Dx) * dzdYSq;
+ const double S3 = (Cy + Cx) * dzdX * dzdY;
+ const double S4 = dzdXSq * dzdXSq;
+ const double S5 = (Dy + Dx) * dzdXSq;
+ const double S6 = Cx * Cy;
+ const double S7 = dydXSq * dydXSq;
+ const double S8 = Dy*dydXSq;
+ const double S9 = dxdXSq * dxdXSq;
+ const double S10 = Dx*dxdXSq;
+ const double S11 = Dy*Dy;
+ const double S12 = Dx*Dy;
+ const double S13 = Dx*Dx;
+ const double S = 2.*(S2+S5+S12)+ 4.*(S7-S8+S9-S10) + 8.*(S3+S6) + S1+S4+S11+S13;
+ const double N = dxdXSq + dxdYSq + dydXSq + dydYSq + dzdXSq + dzdYSq;
+ const double sqrtS = sqrt(S);
+ const double sigma1Sq = 0.5 * (N + sqrtS), sigma2Sq = 0.5 * (N - sqrtS);
+ if (sigma2Sq < Eps) return 0.;
+ const double iCN = 2. / sqrt((sigma1Sq + sigma2Sq) * (1./sigma1Sq + 1./sigma2Sq));
+ if (sign) {
+ const double lnx = dydX*dzdY - dzdX*dydY, lny = dzdX*dxdY - dxdX*dzdY, // Local normal from mapping gradients
+ lnz = dxdX*dydY - dydX*dxdY;
+ const double dp = lnx*nx + lny*ny + lnz*nz; // Dot product to determine element validity
+ if (dp == 0.) return 0.;
+ return (dp > 0.) ? iCN : -iCN;
+ }
+ else
+ return iCN;
+// return std::min(sqrt(sigma1Sq), sqrt(sigma2Sq)) / std::max(sqrt(sigma1Sq), sqrt(sigma2Sq));
+}
+
+
+// Compute the squared Frobenius norm of the inverse of a matrix
+template<bool sign>
+inline double calcInvCondNum3D(double J11, double J12, double J13,
+ double J21, double J22, double J23,
+ double J31, double J32, double J33)
+{
+ const double D = calcDet3x3(J11, J12, J13, J21, J22, J23, J31, J32, J33);
+ if (D == 0.) return 0.;
+ const double I11 = J22*J33-J23*J32, I12 = J13*J32-J12*J33, I13 = J12*J23-J13*J22,
+ I21 = J23*J31-J21*J33, I22 = J11*J33-J13*J31, I23 = J13*J21-J11*J23,
+ I31 = J21*J32-J22*J31, I32 = J12*J31-J11*J32, I33 = J11*J22-J12*J21;
+ const double nSqJ = J11*J11 + J12*J12 + J13*J13
+ + J21*J21 + J22*J22 + J23*J23
+ + J31*J31 + J32*J32 + J33*J33;
+ const double nSqDInvJ = I11*I11 + I12*I12 + I13*I13
+ + I21*I21 + I22*I22 + I23*I23
+ + I31*I31 + I32*I32 + I33*I33;
+ if (sign)
+ return 3.*D/sqrt(nSqJ*nSqDInvJ);
+ else
+ return 3.*std::fabs(D)/sqrt(nSqJ*nSqDInvJ);
+}
+
+
+// Compute condition number and its gradients
+// w.r.t. node positions, at one location in a 2D element
+template<bool sign>
+inline void calcGradInvCondNum2D(double dxdX, double dxdY,
+ double dydX, double dydY,
+ double dzdX, double dzdY,
+ double nx, double ny, double nz,
+ int i, int numMapNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ fullMatrix<double> &IDI)
+{
+ static const double Eps = 1e-10;
+
+ bool posJac = true;
+ if (sign) {
+ const double lnx = dydX*dzdY - dzdX*dydY, lny = dzdX*dxdY - dxdX*dzdY, // Local normal from mapping gradients
+ lnz = dxdX*dydY - dydX*dxdY;
+ const double dp = lnx*nx + lny*ny + lnz*nz; // Dot product to determine element validity
+ posJac = (dp > 0.);
+ }
+
+ const double dxdXSq = dxdX*dxdX, dydXSq = dydX*dydX, dzdXSq = dzdX*dzdX;
+ const double dxdYSq = dxdY*dxdY, dydYSq = dydY*dydY, dzdYSq = dzdY*dzdY;
+ const double Dx = dxdXSq - dxdYSq, Dy = dydXSq - dydYSq;
+ const double Cx = dxdX * dxdY, Cy = dydX * dydY;
+ const double S1 = dzdYSq * dzdYSq;
+ const double S2 = (dzdXSq - Dy - Dx) * dzdYSq;
+ const double S3 = (Cy + Cx) * dzdX * dzdY;
+ const double S4 = dzdXSq * dzdXSq;
+ const double S5 = (Dy + Dx) * dzdXSq;
+ const double S6 = Cx * Cy;
+ const double S7 = dydXSq * dydXSq;
+ const double S8 = Dy*dydXSq;
+ const double S9 = dxdXSq * dxdXSq;
+ const double S10 = Dx*dxdXSq;
+ const double S11 = Dy*Dy;
+ const double S12 = Dx*Dy;
+ const double S13 = Dx*Dx;
+ const double S = 2.*(S2+S5+S12)+ 4.*(S7-S8+S9-S10) + 8.*(S3+S6) + S1+S4+S11+S13;
+ if (S < Eps) { // S == 0. -> Ideal element
+ for (int j = 0; j<3*numMapNodes; j++) IDI(i, j) = 0.;
+ IDI(i, 3*numMapNodes) = posJac ? 1. : -1.;
+ return;
+ }
+ const double N = dxdXSq + dxdYSq + dydXSq + dydYSq + dzdXSq + dzdYSq;
+ const double sqrtS = sqrt(S), invSqrtS = 1./sqrtS;
+ const double sigma1Sq = 0.5 * (N + sqrtS), sigma2Sq = 0.5 * (N - sqrtS);
+// if (sigma2Sq < Eps) { // TODO: Degenerate element (sigma2 == 0.)
+// const double dnx = dydX*nz - dzdX*ny, dny = dzdX*nx - dxdX*nz,
+// dnz = dxdX*ny - dydX*nx;
+//
+// }
+ const double invSigma1Sq = 1./sigma1Sq;
+ const double invSigma2Sq = 1./sigma2Sq;
+ const double invProd = 1. / ((sigma1Sq + sigma2Sq) * (invSigma1Sq + invSigma2Sq));
+ const double invSqrtProd = sqrt(invProd);
+ IDI(i, 3*numMapNodes) = posJac ? 2. * invSqrtProd : -2. * invSqrtProd;
+// std::cout << "DBGTT: Value " << i << "\n";
+// std::cout << "DBGTT: -> I = " << IDI(i, 3*numMapNodes) << "\n";
+// std::cout << "DBGTT: -> S = " << S << ", sqrt(S) = " << sqrtS << "\n";
+// std::cout << "DBGTT: -> N = " << N << "\n";
+// std::cout << "DBGTT: -> sigma1Sq, sigma2Sq = (" << sigma1Sq << ", " << sigma2Sq << ")\n";
+
+ for (int j = 0; j < numMapNodes; j++) {
+ const double &dPhidX = gSMatX(i, j);
+ const double &dPhidY = gSMatY(i, j);
+// std::cout << "DBGTT: -> Node " << j << ":\n";
+
+ const double ddxdXSqdxj = 2.*dPhidX*dxdX, ddxdYSqdxj = 2.*dPhidY*dxdY;
+ const double dDxdxj = ddxdXSqdxj - ddxdYSqdxj;
+ const double dCxdxj = dPhidX*dxdY + dxdX*dPhidY;
+ const double dS2dxj = -dDxdxj * dzdYSq;
+ const double dS3dxj = dCxdxj * dzdX * dzdY;
+ const double dS5dxj = dDxdxj * dzdXSq;
+ const double dS6dxj = dCxdxj * Cy;
+ const double dS9dxj = 2. * ddxdXSqdxj * dxdXSq;
+ const double dS10dxj = dDxdxj*dxdXSq + Dx*ddxdXSqdxj;
+ const double dS12dxj = dDxdxj * Dy;
+ const double dS13dxj = 2. * dDxdxj* Dx;
+ const double dSdxj = 2.*(dS2dxj+dS5dxj+dS12dxj)+ 4.*(dS9dxj-dS10dxj)
+ + 8.*(dS3dxj+dS6dxj) + dS13dxj;
+ const double dNdxj = ddxdXSqdxj + ddxdYSqdxj;
+ const double dsqrtSdxj = 0.5*dSdxj*invSqrtS;
+ const double dsigma1Sqdxj = 0.5 * (dNdxj + dsqrtSdxj),
+ dsigma2Sqdxj = 0.5 * (dNdxj - dsqrtSdxj);
+ const double dinvSigma1Sqdxj = -dsigma1Sqdxj*invSigma1Sq*invSigma1Sq,
+ dinvSigma2Sqdxj = -dsigma2Sqdxj*invSigma2Sq*invSigma2Sq;
+ const double dProddxj = (dsigma1Sqdxj + dsigma2Sqdxj) * (invSigma1Sq + invSigma2Sq)
+ + (sigma1Sq + sigma2Sq) * (dinvSigma1Sqdxj + dinvSigma2Sqdxj);
+ const double mdiCNdxj = dProddxj*invProd*invSqrtProd;
+ IDI(i, j) = posJac ? -mdiCNdxj : mdiCNdxj;
+// std::cout << "DBGTT: -> d(dxdX)dxj = " << dPhidX << ", d(dxdY)dxj = " << dPhidY << "\n";
+// std::cout << "DBGTT: -> dNdxj = " << dNdxj << ", dSdxj = " << dSdxj << "\n";
+// std::cout << "DBGTT: -> dsigma1Sqdxj = " << dsigma1Sqdxj << ", dsigma2Sqdxj = " << dsigma2Sqdxj << "\n";
+// std::cout << "DBGTT: -> dinvSigma1Sqdxj = " << dinvSigma1Sqdxj << ", dinvSigma2Sqdxj = " << dinvSigma2Sqdxj << "\n";
+// std::cout << "DBGTT: -> mdiCNdxj = " << mdiCNdxj
+// << ", -dsigma2Sqdxj/sqrt(sigma1Sq*sigma2Sq) = " << -dsigma2Sqdxj/sqrt(sigma1Sq*sigma2Sq)
+// << ", T1 = " << -dNdxj/sqrt(N*N-S) << ", T21 = " << (0.5*dSdxj/sqrtS)/sqrt(N*N-S) << "\n";
+
+ const double ddydXSqdyj = 2.*dPhidX*dydX, ddydYSqdyj = 2.*dPhidY*dydY;
+ const double dDydyj = ddydXSqdyj - ddydYSqdyj;
+ const double dCydyj = dPhidX*dydY + dydX*dPhidY;
+ const double dS2dyj = -dDydyj * dzdYSq;
+ const double dS3dyj = dCydyj * dzdX * dzdY;
+ const double dS5dyj = dDydyj * dzdXSq;
+ const double dS6dyj = Cx * dCydyj;
+ const double dS7dyj = 2. * ddydXSqdyj * dydXSq;
+ const double dS8dyj = dDydyj*dydXSq + Dy*ddydXSqdyj;
+ const double dS11dyj = 2. * dDydyj * Dy;
+ const double dS12dyj = Dx* dDydyj;
+ const double dSdyj = 2.*(dS2dyj+dS5dyj+dS12dyj)+ 4.*(dS7dyj-dS8dyj)
+ + 8.*(dS3dyj+dS6dyj) + dS11dyj;
+ const double dNdyj = ddydXSqdyj + ddydYSqdyj;
+ const double dsqrtSdyj = 0.5*dSdyj*invSqrtS;
+ const double dsigma1Sqdyj = 0.5 * (dNdyj + dsqrtSdyj),
+ dsigma2Sqdyj = 0.5 * (dNdyj - dsqrtSdyj);
+ const double dinvSigma1Sqdyj = -dsigma1Sqdyj*invSigma1Sq*invSigma1Sq,
+ dinvSigma2Sqdyj = -dsigma2Sqdyj*invSigma2Sq*invSigma2Sq;
+ const double dProddyj = (dsigma1Sqdyj + dsigma2Sqdyj) * (invSigma1Sq + invSigma2Sq)
+ + (sigma1Sq + sigma2Sq) * (dinvSigma1Sqdyj + dinvSigma2Sqdyj);
+ const double mdiCNdyj = dProddyj*invProd*invSqrtProd;
+ IDI(i, j+numMapNodes) = posJac ? -mdiCNdyj : mdiCNdyj;
+// std::cout << "DBGTT: -> d(dydX)dyj = " << dPhidX << ", d(dydY)dyj = " << dPhidY << "\n";
+// std::cout << "DBGTT: -> dNdyj = " << dNdyj << ", dSdyj = " << dSdyj << "\n";
+// std::cout << "DBGTT: -> dsigma1Sqdyj = " << dsigma1Sqdyj << ", dsigma2Sqdyj = " << dsigma2Sqdyj << "\n";
+// std::cout << "DBGTT: -> dinvSigma1Sqdyj = " << dinvSigma1Sqdyj << ", dinvSigma2Sqdyj = " << dinvSigma2Sqdyj << "\n";
+
+ const double ddzdXSqdzj = 2.*dPhidX*dzdX, ddzdYSqdzj = 2.*dPhidY*dzdY;
+ const double dS1dzj = 2. * ddzdYSqdzj * dzdYSq;
+ const double dS2dzj = (dzdXSq - Dy - Dx) * ddzdYSqdzj;
+ const double dS3dzj = (Cy + Cx) *
+ (ddzdXSqdzj*dzdY + dzdX*ddzdYSqdzj);
+ const double dS4dzj = 2. * ddzdXSqdzj * dzdXSq;
+ const double dS5dzj = (Dy + Dx) * ddzdXSqdzj;
+ const double dSdzj = 2.*(dS2dzj+dS5dzj) + 8.*dS3dzj + dS1dzj+dS4dzj;
+ const double dNdzj = ddzdXSqdzj + ddzdYSqdzj;
+ const double dsqrtSdzj = 0.5*dSdzj*invSqrtS;
+ const double dsigma1Sqdzj = 0.5 * (dNdzj + dsqrtSdzj),
+ dsigma2Sqdzj = 0.5 * (dNdzj - dsqrtSdzj);
+ const double dinvSigma1Sqdzj = -dsigma1Sqdzj*invSigma1Sq*invSigma1Sq,
+ dinvSigma2Sqdzj = -dsigma2Sqdzj*invSigma2Sq*invSigma2Sq;
+ const double dProddzj = (dsigma1Sqdzj + dsigma2Sqdzj) * (invSigma1Sq + invSigma2Sq)
+ + (sigma1Sq + sigma2Sq) * (dinvSigma1Sqdzj + dinvSigma2Sqdzj);
+ const double mdiCNdzj = dProddzj*invProd*invSqrtProd;
+ IDI(i, j+2*numMapNodes) = posJac ? -mdiCNdzj : mdiCNdzj;
+// std::cout << "DBGTT: -> d(dzdX)dzj = " << dPhidX << ", d(dzdY)dzj = " << dPhidY << "\n";
+// std::cout << "DBGTT: -> dNdzj = " << dNdzj << ", dSdzj = " << dSdzj << "\n";
+// std::cout << "DBGTT: -> dsigma1Sqdzj = " << dsigma1Sqdzj << ", dsigma2Sqdzj = " << dsigma2Sqdzj << "\n";
+// std::cout << "DBGTT: -> dinvSigma1Sqdzj = " << dinvSigma1Sqdzj << ", dinvSigma2Sqdzj = " << dinvSigma2Sqdzj << "\n";
+// std::cout << "DBGTT: -> dId{x,y,z} = (" << IDI(i, j) << ", "
+// << IDI(i, j+numMapNodes) << ", " << IDI(i, j+2*numMapNodes) << ")\n";
+ }
+}
+
+
+// Compute condition number and its gradients
+// w.r.t. node positions, at one location in a 3D element
+template<bool sign>
+inline void calcGradInvCondNum3D(double dxdX, double dxdY, double dxdZ,
+ double dydX, double dydY, double dydZ,
+ double dzdX, double dzdY, double dzdZ,
+ int i, int numMapNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ fullMatrix<double> &IDI)
+{
+ const double normJSq = dxdX*dxdX + dxdY*dxdY + dxdZ*dxdZ
+ + dydX*dydX + dydY*dydY + dydZ*dydZ
+ + dzdX*dzdX + dzdY*dzdY + dzdZ*dzdZ;
+ const double I11 = dydY*dzdZ - dydZ*dzdY, I12 = dxdZ*dzdY - dxdY*dzdZ,
+ I13 = dxdY*dydZ - dxdZ*dydY, I21 = dydZ*dzdX - dydX*dzdZ,
+ I22 = dxdX*dzdZ - dxdZ*dzdX, I23 = dxdZ*dydX - dxdX*dydZ,
+ I31 = dydX*dzdY - dydY*dzdX, I32 = dxdY*dzdX - dxdX*dzdY,
+ I33 = dxdX*dydY - dxdY*dydX;
+ const double normISq = I11*I11 + I12*I12 + I13*I13
+ + I21*I21 + I22*I22 + I23*I23
+ + I31*I31 + I32*I32 + I33*I33;
+ const double invProd = 1./(normJSq*normISq), invSqrtProd = sqrt(invProd);
+ const double D = calcDet3x3(dxdX, dxdY, dxdZ,
+ dydX, dydY, dydZ,
+ dzdX, dzdY, dzdZ);
+ const bool reverse = (!sign && (D < 0.));
+ const double sICN = 3.*D*invSqrtProd;
+ IDI(i, 3*numMapNodes) = reverse ? -sICN : sICN;
+
+ for (int j = 0; j < numMapNodes; j++) {
+ const double &dPhidX = gSMatX(i, j);
+ const double &dPhidY = gSMatY(i, j);
+ const double &dPhidZ = gSMatZ(i, j);
+
+ const double dNormJSqdxj = 2.*(dPhidX*dxdX + dPhidY*dxdY + dPhidZ*dxdZ);
+ const double dNormISqdxj = 2.*((dPhidZ*dzdY - dPhidY*dzdZ)*I12 + (dPhidY*dydZ - dPhidZ*dydY)*I13
+ + (dPhidX*dzdZ - dPhidZ*dzdX)*I22 + (dPhidZ*dydX - dPhidX*dydZ)*I23
+ + (dPhidY*dzdX - dPhidX*dzdY)*I32 + (dPhidX*dydY - dPhidY*dydX)*I33);
+ const double dProddxj = dNormJSqdxj*normISq + dNormISqdxj*normJSq;
+ const double dDdxj = dPhidX * dydY * dzdZ + dzdX * dPhidY * dydZ +
+ dydX * dzdY * dPhidZ - dzdX * dydY * dPhidZ -
+ dPhidX * dzdY * dydZ - dydX * dPhidY * dzdZ;
+ const double dsICNdxj = 3. * (dDdxj*invSqrtProd - 0.5*D*dProddxj*invProd*invSqrtProd);
+ IDI(i, j) = reverse ? -dsICNdxj : dsICNdxj;
+
+ const double dNormJSqdyj = 2.*(dPhidX*dydX + dPhidY*dydY + dPhidZ*dydZ);
+ const double dNormISqdyj = 2.*((dPhidY*dzdZ - dPhidZ*dzdY)*I11 + (dxdY*dPhidZ - dxdZ*dPhidY)*I13
+ + (dPhidZ*dzdX - dPhidX*dzdZ)*I21 + (dxdZ*dPhidX - dxdX*dPhidZ)*I23
+ + (dPhidX*dzdY - dPhidY*dzdX)*I31 + (dxdX*dPhidY - dxdY*dPhidX)*I33);
+ const double dProddyj = dNormJSqdyj*normISq + dNormISqdyj*normJSq;
+ const double dDdyj = dxdX * dPhidY * dzdZ + dzdX * dxdY * dPhidZ +
+ dPhidX * dzdY * dxdZ - dzdX * dPhidY * dxdZ -
+ dxdX * dzdY * dPhidZ - dPhidX * dxdY * dzdZ;
+ const double dsICNdyj = 3. * (dDdyj*invSqrtProd - 0.5*D*dProddyj*invProd*invSqrtProd);
+ IDI(i, j+numMapNodes) = reverse ? -dsICNdyj : dsICNdyj;
+
+ const double dNormJSqdzj = 2.*(dPhidX*dzdX + dPhidY*dzdY + dPhidZ*dzdZ);
+ const double dNormISqdzj = 2.*((dydY*dPhidZ - dydZ*dPhidY)*I11 + (dxdZ*dPhidY - dxdY*dPhidZ)*I12
+ + (dydZ*dPhidX - dydX*dPhidZ)*I21 + (dxdX*dPhidZ - dxdZ*dPhidX)*I22
+ + (dydX*dPhidY - dydY*dPhidX)*I31 + (dxdY*dPhidX - dxdX*dPhidY)*I32);
+ const double dProddzj = dNormJSqdzj*normISq + dNormISqdzj*normJSq;
+ const double dDdzj = dxdX * dydY * dPhidZ + dPhidX * dxdY * dydZ +
+ dydX * dPhidY * dxdZ - dPhidX * dydY * dxdZ -
+ dxdX * dPhidY * dydZ - dydX * dxdY * dPhidZ;
+ const double dsICNdzj = 3. * (dDdzj*invSqrtProd - 0.5*D*dProddzj*invProd*invSqrtProd);
+ IDI(i, j+2*numMapNodes) = reverse ? -dsICNdzj : dsICNdzj;
+ }
+}
+
+
+}
+
+
+CondNumBasis::CondNumBasis(int tag, int cnOrder) :
+ _tag(tag), _dim(ElementType::DimensionFromTag(tag)),
+ _condNumOrder(cnOrder >= 0 ? cnOrder : condNumOrder(tag))
+{
+ const bool serendip = false;
+ FuncSpaceData data(true, tag, _condNumOrder, &serendip);
+
+ _gradBasis = BasisFactory::getGradientBasis(data);
+
+ fullMatrix<double> lagPoints; // Sampling points
+ gmshGeneratePoints(data, lagPoints);
+ _nCondNumNodes = lagPoints.size1();
+ _nMapNodes = BasisFactory::getNodalBasis(tag)->getNumShapeFunctions();
+
+ // Store shape function gradients of mapping at condition number nodes
+ _gradBasis = BasisFactory::getGradientBasis(tag, _condNumOrder);
+
+ // Compute shape function gradients of primary mapping at barycenter,
+ // in order to compute normal to straight element
+ const int parentType = ElementType::ParentTypeFromTag(tag);
+ const int primMapType = ElementType::getTag(parentType, 1, false);
+ const nodalBasis *primMapBasis = BasisFactory::getNodalBasis(primMapType);
+ _nPrimMapNodes = primMapBasis->getNumShapeFunctions();
+
+ double xBar = 0., yBar = 0., zBar = 0.;
+ double barycenter[3] = {0., 0., 0.};
+ for (int i = 0; i < _nPrimMapNodes; i++) {
+ for (int j = 0; j < primMapBasis->points.size2(); ++j) {
+ barycenter[j] += primMapBasis->points(i, j);
+ }
+ }
+ barycenter[0] /= _nPrimMapNodes;
+ barycenter[1] /= _nPrimMapNodes;
+ barycenter[2] /= _nPrimMapNodes;
+
+ double (*barDPsi)[3] = new double[_nPrimMapNodes][3];
+ primMapBasis->df(xBar, yBar, zBar, barDPsi);
+
+ // TODO: Make primGradShape from ideal element
+ primGradShapeBarycenterX.resize(_nPrimMapNodes);
+ primGradShapeBarycenterY.resize(_nPrimMapNodes);
+ primGradShapeBarycenterZ.resize(_nPrimMapNodes);
+ for (int j=0; j<_nPrimMapNodes; j++) {
+ primGradShapeBarycenterX(j) = barDPsi[j][0];
+ primGradShapeBarycenterY(j) = barDPsi[j][1];
+ primGradShapeBarycenterZ(j) = barDPsi[j][2];
+ }
+
+ delete[] barDPsi;
+}
+
+
+int CondNumBasis::condNumOrder(int tag)
+{
+ const int parentType = ElementType::ParentTypeFromTag(tag);
+ const int order = ElementType::OrderFromTag(tag);
+ return condNumOrder(parentType, order);
+}
+
+
+int CondNumBasis::condNumOrder(int parentType, int order)
+{
+ switch (parentType) {
+ case TYPE_PNT : return 0;
+ case TYPE_LIN : return order - 1;
+ case TYPE_TRI : return 2*order - 2;
+ case TYPE_QUA : return 2*order - 1;
+ case TYPE_TET : return 3*order - 3;
+ case TYPE_PRI : return 3*order - 1;
+ case TYPE_HEX : return 3*order - 1;
+ case TYPE_PYR : return 3*order - 3;
+ default :
+ Msg::Error("Unknown element type %d, return order 0", parentType);
+ return 0;
+ }
+}
+
+
+// Calculate the inverse condition number in Frobenius norm for one element, with normal vectors to straight
+// element for regularization. Evaluation points depend on the given matrices for shape function gradients.
+template<bool sign>
+inline void CondNumBasis::getInvCondNumGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullVector<double> &condNum) const
+{
+ switch (_dim) {
+
+ case 0 : {
+ for (int i = 0; i < nCondNumNodes; i++) condNum(i) = 1.;
+ break;
+ }
+
+ case 1 : {
+ Msg::Warning("Inverse condition number not implemented in 1D");
+ condNum.setAll(0.);
+ break;
+ }
+
+ case 2 : {
+ fullMatrix<double> dxyzdX(nCondNumNodes,3), dxyzdY(nCondNumNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ for (int i = 0; i < nCondNumNodes; 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 &nx = normals(0,0), &ny = normals(0,1), &nz = normals(0,2);
+ condNum(i) = calcInvCondNum2D<sign>(dxdX, dxdY, dydX, dydY, dzdX, dzdY, nx, ny, nz);
+ }
+ break;
+ }
+
+ case 3 : {
+ fullMatrix<double> dxyzdX(nCondNumNodes, 3), dxyzdY(nCondNumNodes, 3), dxyzdZ(nCondNumNodes, 3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ gSMatZ.mult(nodesXYZ, dxyzdZ);
+ for (int i = 0; i < nCondNumNodes; 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);
+ condNum(i) = calcInvCondNum3D<sign>(dxdX, dxdY, dxdZ,
+ dydX, dydY, dydZ,
+ dzdX, dzdY, dzdZ);
+ }
+ break;
+ }
+ }
+}
+
+
+void CondNumBasis::getInvCondNumGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &invCond) const
+{
+ fullMatrix<double> dumNormals;
+ getInvCondNumGeneral<false>(nCondNumNodes, gSMatX, gSMatY, gSMatZ,
+ nodesXYZ, dumNormals, invCond);
+}
+
+
+void CondNumBasis::getSignedInvCondNumGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullVector<double> &invCond) const
+{
+ getInvCondNumGeneral<true>(nCondNumNodes, gSMatX, gSMatY, gSMatZ,
+ nodesXYZ, normals, invCond);
+}
+
+
+// Calculate the inverse condition number in Frobenius norm and its gradients w.r.t. node position,
+// with normal vectors to straight element for regularization. Evaluation points depend on the
+// given matrices for shape function gradients.
+template<bool sign>
+inline void CondNumBasis::getInvCondNumAndGradientsGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &IDI) const
+{
+ fullMatrix<double> JDJ(nCondNumNodes, 3*_nMapNodes+1);
+
+ switch (_dim) {
+
+ case 0 : {
+ for (int i = 0; i < nCondNumNodes; i++) {
+ for (int j = 0; j < _nMapNodes; j++) {
+ IDI(i,j) = 0.;
+ IDI(i,j+1*_nMapNodes) = 0.;
+ IDI(i,j+2*_nMapNodes) = 0.;
+ }
+ IDI(i,3*_nMapNodes) = 1.;
+ }
+ break;
+ }
+
+ case 1 : {
+ Msg::Warning("Inverse condition number not implemented in 1D");
+ IDI.setAll(0.);
+ break;
+ }
+
+ case 2 : {
+ fullMatrix<double> dxyzdX(nCondNumNodes,3), dxyzdY(nCondNumNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ for (int i = 0; i < nCondNumNodes; 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 &nx = normals(0,0), &ny = normals(0,1), &nz = normals(0,2);
+ calcGradInvCondNum2D<sign>(dxdX, dxdY, dydX, dydY, dzdX, dzdY,
+ nx, ny, nz, i, _nMapNodes, gSMatX, gSMatY, IDI);
+ }
+ break;
+ }
+
+ case 3 : {
+ fullMatrix<double> dxyzdX(nCondNumNodes,3), dxyzdY(nCondNumNodes,3), dxyzdZ(nCondNumNodes,3);
+ gSMatX.mult(nodesXYZ, dxyzdX);
+ gSMatY.mult(nodesXYZ, dxyzdY);
+ gSMatZ.mult(nodesXYZ, dxyzdZ);
+ for (int i = 0; i < nCondNumNodes; 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);
+ calcGradInvCondNum3D<sign>(dxdX, dxdY, dxdZ,
+ dydX, dydY, dydZ,
+ dzdX, dzdY, dzdZ,
+ i, _nMapNodes,
+ gSMatX, gSMatY, gSMatZ, IDI);
+ }
+ break;
+ }
+ }
+}
+
+
+void CondNumBasis::getInvCondNumAndGradientsGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ fullMatrix<double> &IDI) const
+{
+ fullMatrix<double> dumNormals;
+ getInvCondNumAndGradientsGeneral<false>(nCondNumNodes, gSMatX, gSMatY, gSMatZ,
+ nodesXYZ, dumNormals, IDI);
+}
+
+
+void CondNumBasis::getSignedInvCondNumAndGradientsGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &IDI) const
+{
+ getInvCondNumAndGradientsGeneral<true>(nCondNumNodes, gSMatX, gSMatY, gSMatZ,
+ nodesXYZ, normals, IDI);
+}
diff --git a/Numeric/CondNumBasis.h b/Numeric/CondNumBasis.h
new file mode 100644
index 0000000000000000000000000000000000000000..bcbd98d8ac8351ef9068ac2d7c971d5fddfd7824
--- /dev/null
+++ b/Numeric/CondNumBasis.h
@@ -0,0 +1,125 @@
+// Gmsh - Copyright (C) 1997-2014 C. Geuzaine, J.-F. Remacle
+//
+// See the LICENSE.txt file for license information. Please report all
+// bugs and problems to the public mailing list <gmsh@geuz.org>.
+
+#ifndef _CONDNUMBASIS_H_
+#define _CONDNUMBASIS_H_
+
+#include <map>
+#include <vector>
+#include "JacobianBasis.h"
+#include "fullMatrix.h"
+
+
+class CondNumBasis {
+ private:
+ const GradientBasis *_gradBasis;
+
+ const int _tag, _dim, _condNumOrder;
+
+ fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY,
+ primGradShapeBarycenterZ;
+
+ int _nCondNumNodes;
+ int _nMapNodes, _nPrimMapNodes;
+
+ public :
+ CondNumBasis(int tag, int cnOrder = -1);
+
+ // Get methods
+ inline int getCondNumOrder() const { return _condNumOrder; }
+ inline int getNumCondNumNodes() const { return _nCondNumNodes; }
+ inline int getNumMapNodes() const { return _nMapNodes; }
+ inline int getNumPrimMapNodes() const { return _nPrimMapNodes; }
+
+ // Order calculation
+ static int condNumOrder(int tag);
+ static int condNumOrder(int parentType, int order);
+
+ // Condition number evaluation methods
+ inline void getInvCondNum(const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &invCond) const {
+ getInvCondNumGeneral(_nCondNumNodes,
+ _gradBasis->gradShapeIdealMatX,
+ _gradBasis->gradShapeIdealMatY,
+ _gradBasis->gradShapeIdealMatZ,
+ nodesXYZ, invCond);
+ }
+ inline void getSignedInvCondNum(const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullVector<double> &invCond) const {
+ getSignedInvCondNumGeneral(_nCondNumNodes,
+ _gradBasis->gradShapeIdealMatX,
+ _gradBasis->gradShapeIdealMatY,
+ _gradBasis->gradShapeIdealMatZ,
+ nodesXYZ, normals, invCond);
+ }
+ inline void getInvCondNumAndGradients(const fullMatrix<double> &nodesXYZ,
+ fullMatrix<double> &IDI) const {
+ getInvCondNumAndGradientsGeneral(_nCondNumNodes,
+ _gradBasis->gradShapeIdealMatX,
+ _gradBasis->gradShapeIdealMatY,
+ _gradBasis->gradShapeIdealMatZ,
+ nodesXYZ, IDI);
+ }
+ inline void getSignedInvCondNumAndGradients(const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &IDI) const {
+ getSignedInvCondNumAndGradientsGeneral(_nCondNumNodes,
+ _gradBasis->gradShapeIdealMatX,
+ _gradBasis->gradShapeIdealMatY,
+ _gradBasis->gradShapeIdealMatZ,
+ nodesXYZ, normals, IDI);
+ }
+
+
+ private :
+ template<bool sign>
+ void getInvCondNumGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullVector<double> &invCond) const;
+ void getInvCondNumGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ fullVector<double> &invCond) const;
+ void getSignedInvCondNumGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullVector<double> &invCond) const;
+
+ template<bool sign>
+ void getInvCondNumAndGradientsGeneral(int nCondNumNodes,
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &IDI) const;
+ void getInvCondNumAndGradientsGeneral(int nCondNumNodes, // No normal given -> unsigned measure
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ fullMatrix<double> &IDI) const;
+ void getSignedInvCondNumAndGradientsGeneral(int nCondNumNodes, // Normals given -> signed measure
+ const fullMatrix<double> &gSMatX,
+ const fullMatrix<double> &gSMatY,
+ const fullMatrix<double> &gSMatZ,
+ const fullMatrix<double> &nodesXYZ,
+ const fullMatrix<double> &normals,
+ fullMatrix<double> &IDI) const;
+
+
+};
+
+#endif
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 6a57c6313474276f7d75fc1b5dc476c2aa49376a..d67dcad36f4b4a2e2decb8fd059982fa1d46ebf5 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -22,30 +22,6 @@ inline double calcDet3D(double M11, double M12, double M13,
+ M13 * (M21*M32 - M22*M31);
}
-// Compute the squared Frobenius norm a matrix
-inline double calcFrobNormSq3D(double M11, double M12, double M13,
- double M21, double M22, double M23,
- double M31, double M32, double M33)
-{
- return M11*M11 + M12*M12 + M13*M13
- + M21*M21 + M22*M22 + M23*M23
- + M31*M31 + M32*M32 + M33*M33;
-}
-
-// Compute the squared Frobenius norm of the inverse of a matrix
-inline double calcFrobNormSqInv3D(double M11, double M12, double M13,
- double M21, double M22, double M23,
- double M31, double M32, double M33)
-{
- const double invD = 1./calcDet3D(M11, M12, M13, M21, M22, M23, M31, M32, M33);
- const double I11 = M22*M33-M23*M32, I12 = M13*M32-M12*M33, I13 = M12*M23-M13*M22,
- I21 = M23*M31-M21*M33, I22 = M11*M33-M13*M31, I23 = M13*M21-M11*M23,
- I31 = M21*M32-M22*M31, I32 = M12*M31-M11*M32, I33 = M11*M22-M12*M21;
- return (I11*I11 + I12*I12 + I13*I13
- + I21*I21 + I22*I22 + I23*I23
- + I31*I31 + I32*I32 + I33*I33)*invD*invD;
-}
-
// Compute signed Jacobian and its gradients w.r.t.
// node positions, at one location in a 1D element
inline void calcJDJ1D(double dxdX, double dxdY, double dxdZ,
@@ -129,156 +105,6 @@ inline void calcJDJ3D(double dxdX, double dxdY, double dxdZ,
dzdX, dzdY, dzdZ);
}
-// Compute condition number and its gradients
-// w.r.t. node positions, at one location in a 1D element
-// FIXME: TO BE UPDATED
-inline void calcIDI1D(double dxdX, double dxdY, double dxdZ,
- double dydX, double dydY, double dydZ,
- double dzdX, double dzdY, double dzdZ,
- int i, int numMapNodes,
- const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &JDJ,
- fullMatrix<double> &IDI)
-{
- const double &J = JDJ(i, 3*numMapNodes), JSq = J*J, JQuadInv = 1./(JSq*JSq);
- const double A1 = dydY*dzdZ - dydZ*dzdY, A2 = dydX*dzdZ - dydZ*dzdX,
- A3 = dxdY*dzdZ - dxdZ*dzdY, A4 = dxdX*dzdZ - dxdZ*dzdX,
- A5 = dydX*dzdY - dydY*dzdX, A6 = dxdX*dzdY - dxdY*dzdX,
- A7 = dxdY*dydZ - dxdZ*dydY, A8 = dxdX*dydZ - dxdZ*dydX,
- A9 = dxdX*dydY - dxdY*dydX;
- const double A = A1*A1 + A2*A2 + A3*A3 + A4*A4 + A5*A5 + A6*A6 + A7*A7 + A8*A8 + A9*A9;
- const double nInvJSq = A/JSq;
- const double nJSq = dxdX*dxdX + dxdY*dxdY + dxdZ*dxdZ
- + dydX*dydX + dydY*dydY + dydZ*dydZ
- + dzdX*dzdX + dzdY*dzdY + dzdZ*dzdZ;
- const double invProd = 1./(nJSq*nInvJSq), sqrtInvProd = sqrt(invProd);
- IDI(i, 3*numMapNodes) = 3.*sqrtInvProd;
- for (int j = 0; j < numMapNodes; j++) {
- const double &dPhidX = gSMatX(i, j);
- const double &dJdxj = JDJ(i, j), &dJdyj = JDJ(i, j+numMapNodes),
- &dJdzj = JDJ(i, j+2*numMapNodes);
- const double dAdxj = 2.*(dPhidX*dzdZ * A4 + dPhidX*dzdY * A6
- + dPhidX*dydZ * A8 + dPhidX*dydY * A9);
- const double dnInvJSqdxj = (dAdxj*JSq-2.*dJdxj*J*A)*JQuadInv;
- const double dnJSqdxj = 2.*dPhidX*dxdX;
- IDI(i, j) = -1.5 * (dnJSqdxj*nInvJSq+nJSq*dnInvJSqdxj)*invProd*sqrtInvProd;
- const double dAdyj = 2.*(dPhidX*dzdZ * A2 + dPhidX*dzdY * A5
- - dxdZ*dPhidX * A8 - dxdY*dPhidX * A9);
- const double dnInvJSqdyj = (dAdyj*JSq-2.*dJdyj*J*A)*JQuadInv;
- const double dnJSqdyj = 2.*dPhidX*dydX;
- IDI(i, j+numMapNodes) = -1.5 * (dnJSqdyj*nInvJSq+nJSq*dnInvJSqdyj)*invProd*sqrtInvProd;
- const double dAdzj = 2.*(-dydZ*dPhidX * A2 - dxdZ*dPhidX * A4
- - dydY*dPhidX * A5 - dxdY*dPhidX * A6);
- const double dnInvJSqdzj = (dAdzj*JSq-2.*dJdzj*J*A)*JQuadInv;
- const double dnJSqdzj = 2.*dPhidX*dzdX;
- IDI(i, j+2*numMapNodes) = -1.5 * (dnJSqdzj*nInvJSq + nJSq*dnInvJSqdzj) * invProd*sqrtInvProd;
- }
-}
-
-// Compute condition number and its gradients
-// w.r.t. node positions, at one location in a 2D element
-// FIXME: TO BE UPDATED
-inline void calcIDI2D(double dxdX, double dxdY, double dxdZ,
- double dydX, double dydY, double dydZ,
- double dzdX, double dzdY, double dzdZ,
- int i, int numMapNodes,
- const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY,
- const fullMatrix<double> &JDJ,
- fullMatrix<double> &IDI)
-{
- const double &J = JDJ(i, 3*numMapNodes), JSq = J*J, JQuadInv = 1./(JSq*JSq);
- const double A1 = dydY*dzdZ - dydZ*dzdY, A2 = dydX*dzdZ - dydZ*dzdX,
- A3 = dxdY*dzdZ - dxdZ*dzdY, A4 = dxdX*dzdZ - dxdZ*dzdX,
- A5 = dydX*dzdY - dydY*dzdX, A6 = dxdX*dzdY - dxdY*dzdX,
- A7 = dxdY*dydZ - dxdZ*dydY, A8 = dxdX*dydZ - dxdZ*dydX,
- A9 = dxdX*dydY - dxdY*dydX;
- const double A = A1*A1 + A2*A2 + A3*A3 + A4*A4 + A5*A5 + A6*A6 + A7*A7 + A8*A8 + A9*A9;
- const double nInvJSq = A/JSq;
- const double nJSq = dxdX*dxdX + dxdY*dxdY + dxdZ*dxdZ
- + dydX*dydX + dydY*dydY + dydZ*dydZ
- + dzdX*dzdX + dzdY*dzdY + dzdZ*dzdZ;
- const double invProd = 1./(nJSq*nInvJSq), sqrtInvProd = sqrt(invProd);
- IDI(i, 3*numMapNodes) = 3.*sqrtInvProd;
- for (int j = 0; j < numMapNodes; j++) {
- const double &dPhidX = gSMatX(i, j);
- const double &dPhidY = gSMatY(i, j);
- const double &dJdxj = JDJ(i, j), &dJdyj = JDJ(i, j+numMapNodes),
- &dJdzj = JDJ(i, j+2*numMapNodes);
- const double dAdxj = 2.*(dPhidY*dzdZ * A3 + dPhidX*dzdZ * A4
- + (dPhidX*dzdY - dPhidY*dzdX) * A6 + dPhidY*dydZ * A7
- + dPhidX*dydZ * A8 + (dPhidX*dydY - dPhidY*dydX) * A9);
- const double dnInvJSqdxj = (dAdxj*JSq-2.*dJdxj*J*A)*JQuadInv;
- const double dnJSqdxj = 2.*(dPhidX*dxdX + dPhidY*dxdY);
- IDI(i, j) = -1.5 * (dnJSqdxj*nInvJSq+nJSq*dnInvJSqdxj)*invProd*sqrtInvProd;
- const double dAdyj = 2.*(dPhidY*dzdZ * A1 + dPhidX*dzdZ * A2
- + (dPhidX*dzdY - dPhidY*dzdX) * A5 - dxdZ*dPhidY * A7
- - dxdZ*dPhidX * A8 + (dxdX*dPhidY - dxdY*dPhidX) * A9);
- const double dnInvJSqdyj = (dAdyj*JSq-2.*dJdyj*J*A)*JQuadInv;
- const double dnJSqdyj = 2.*(dPhidX*dydX + dPhidY*dydY);
- IDI(i, j+numMapNodes) = -1.5 * (dnJSqdyj*nInvJSq+nJSq*dnInvJSqdyj)*invProd*sqrtInvProd;
- const double dAdzj = 2.*(-dydZ*dPhidY * A1 - dydZ*dPhidX * A2
- - dxdZ*dPhidY * A3 - dxdZ*dPhidX * A4
- + (dydX*dPhidY - dydY*dPhidX) * A5 + (dxdX*dPhidY - dxdY*dPhidX) * A6);
- const double dnInvJSqdzj = (dAdzj*JSq-2.*dJdzj*J*A)*JQuadInv;
- const double dnJSqdzj = 2.*(dPhidX*dzdX + dPhidY*dzdY);
- IDI(i, j+2*numMapNodes) = -1.5 * (dnJSqdzj*nInvJSq + nJSq*dnInvJSqdzj) * invProd*sqrtInvProd;
- }
-}
-
-// Compute condition number and its gradients
-// w.r.t. node positions, at one location in a 3D element
-// FIXME: TO BE UPDATED
-inline void calcIDI3D(double dxdX, double dxdY, double dxdZ,
- double dydX, double dydY, double dydZ,
- double dzdX, double dzdY, double dzdZ,
- int i, int numMapNodes,
- const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY,
- const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &JDJ,
- fullMatrix<double> &IDI)
-{
- const double &J = JDJ(i, 3*numMapNodes), JSq = J*J, JQuadInv = 1./(JSq*JSq);
- const double A1 = dydY*dzdZ - dydZ*dzdY, A2 = dydX*dzdZ - dydZ*dzdX,
- A3 = dxdY*dzdZ - dxdZ*dzdY, A4 = dxdX*dzdZ - dxdZ*dzdX,
- A5 = dydX*dzdY - dydY*dzdX, A6 = dxdX*dzdY - dxdY*dzdX,
- A7 = dxdY*dydZ - dxdZ*dydY, A8 = dxdX*dydZ - dxdZ*dydX,
- A9 = dxdX*dydY - dxdY*dydX;
- const double A = A1*A1 + A2*A2 + A3*A3 + A4*A4 + A5*A5 + A6*A6 + A7*A7 + A8*A8 + A9*A9;
- const double nInvJSq = A/JSq;
- const double nJSq = dxdX*dxdX + dxdY*dxdY + dxdZ*dxdZ
- + dydX*dydX + dydY*dydY + dydZ*dydZ
- + dzdX*dzdX + dzdY*dzdY + dzdZ*dzdZ;
- const double invProd = 1./(nJSq*nInvJSq), sqrtInvProd = sqrt(invProd);
- IDI(i, 3*numMapNodes) = 3.*sqrtInvProd;
- for (int j = 0; j < numMapNodes; j++) {
- const double &dPhidX = gSMatX(i, j);
- const double &dPhidY = gSMatY(i, j);
- const double &dPhidZ = gSMatZ(i, j);
- const double &dJdxj = JDJ(i, j), &dJdyj = JDJ(i, j+numMapNodes),
- &dJdzj = JDJ(i, j+2*numMapNodes);
- const double dAdxj = 2.*((dPhidY*dzdZ - dPhidZ*dzdY) * A3 + (dPhidX*dzdZ - dPhidZ*dzdX) * A4
- + (dPhidX*dzdY - dPhidY*dzdX) * A6 + (dPhidY*dydZ - dPhidZ*dydY) * A7
- + (dPhidX*dydZ - dPhidZ*dydX) * A8 + (dPhidX*dydY - dPhidY*dydX) * A9);
- const double dnInvJSqdxj = (dAdxj*JSq-2.*dJdxj*J*A)*JQuadInv;
- const double dnJSqdxj = 2.*(dPhidX*dxdX + dPhidY*dxdY + dPhidZ*dxdZ);
- IDI(i, j) = -1.5*(dnJSqdxj*nInvJSq+nJSq*dnInvJSqdxj)*invProd*sqrtInvProd;
- const double dAdyj = 2.*((dPhidY*dzdZ - dPhidZ*dzdY) * A1 + (dPhidX*dzdZ - dPhidZ*dzdX) * A2
- + (dPhidX*dzdY - dPhidY*dzdX) * A5 + (dxdY*dPhidZ - dxdZ*dPhidY) * A7
- + (dxdX*dPhidZ - dxdZ*dPhidX) * A8 + (dxdX*dPhidY - dxdY*dPhidX) * A9);
- const double dnInvJSqdyj = (dAdyj*JSq-2.*dJdyj*J*A)*JQuadInv;
- const double dnJSqdyj = 2.*(dPhidX*dydX + dPhidY*dydY + dPhidZ*dydZ);
- IDI(i, j+numMapNodes) = -1.5*(dnJSqdyj*nInvJSq+nJSq*dnInvJSqdyj)*invProd*sqrtInvProd;
- const double dAdzj = 2.*((dydY*dPhidZ - dydZ*dPhidY) * A1 + (dydX*dPhidZ - dydZ*dPhidX) * A2
- + (dxdY*dPhidZ - dxdZ*dPhidY) * A3 + (dxdX*dPhidZ - dxdZ*dPhidX) * A4
- + (dydX*dPhidY - dydY*dPhidX) * A5 + (dxdX*dPhidY - dxdY*dPhidX) * A6);
- const double dnInvJSqdzj = (dAdzj*JSq-2.*dJdzj*J*A)*JQuadInv;
- const double dnJSqdzj = 2.*(dPhidX*dzdX + dPhidY*dzdY + dPhidZ*dzdZ);
- IDI(i, j+2*numMapNodes) = -1.5 * (dnJSqdzj*nInvJSq + nJSq*dnInvJSqdzj) * invProd*sqrtInvProd;
- }
-}
-
}
GradientBasis::GradientBasis(FuncSpaceData data)
@@ -304,6 +130,12 @@ GradientBasis::GradientBasis(FuncSpaceData data)
gradShapeMatZ(i, j) = allDPsi(j, 3*i+2);
}
}
+
+ gradShapeIdealMatX = gradShapeMatX;
+ gradShapeIdealMatY = gradShapeMatY;
+ gradShapeIdealMatZ = gradShapeMatZ;
+ mapFromIdealElement(_data.elementType(), &gradShapeIdealMatX,
+ &gradShapeIdealMatY, &gradShapeIdealMatZ);
}
void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
@@ -316,6 +148,16 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
}
+void GradientBasis::getIdealGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> *dxyzdX,
+ fullMatrix<double> *dxyzdY,
+ fullMatrix<double> *dxyzdZ) const
+{
+ if (dxyzdX) gradShapeIdealMatX.mult(nodes, *dxyzdX);
+ if (dxyzdY) gradShapeIdealMatY.mult(nodes, *dxyzdY);
+ if (dxyzdZ) gradShapeIdealMatZ.mult(nodes, *dxyzdZ);
+}
+
void GradientBasis::mapFromIdealElement(int type,
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
@@ -358,6 +200,7 @@ void GradientBasis::mapFromIdealElement(int type,
dxyzdZ->scale(cTet[2]);
dxyzdZ->axpy(*dxyzdX, cTet[0]);
dxyzdZ->axpy(*dxyzdY, cTet[1]);
+ break;
}
}
}
@@ -868,211 +711,6 @@ inline void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
}
-// Calculate the inverse condition number in Frobenius norm 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::getCondNumGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &invCond) const
-{
- switch (_dim) {
-
- case 0 : {
- for (int i = 0; i < nJacNodes; i++) invCond(i) = 1.;
- break;
- }
-
- case 1 : {
- fullMatrix<double> normals(2,3);
- getPrimNormals1D(nodesXYZ,normals);
- fullMatrix<double> dxyzdX(nJacNodes,3);
- gSMatX.mult(nodesXYZ, dxyzdX);
- for (int i = 0; i < nJacNodes; i++) {
- const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
- const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
- const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
- const double nJSq = calcFrobNormSq3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ);
- const double nInvJSq = calcFrobNormSqInv3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ);
- invCond(i) = 3./sqrt(nJSq*nInvJSq);
- }
- break;
- }
-
- case 2 : {
- fullMatrix<double> normal(1,3);
- const double scaleSq = getPrimNormal2D(nodesXYZ, normal, ideal), scale = sqrt(scaleSq);
- normal(0, 0) *= scale; normal(0, 1) *= scale; normal(0, 2) *= scale;
- 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);
- const double &dxdZ = normal(0,0), &dydZ = normal(0,1), &dzdZ = normal(0,2);
- const double nJSq = calcFrobNormSq3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ);
- const double nInvJSq = calcFrobNormSqInv3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ);
- invCond(i) = 3./sqrt(nJSq*nInvJSq);
- }
- break;
- }
-
- case 3 : {
- fullMatrix<double> dum;
- fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3), dxyzdZ(nJacNodes,3);
- 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);
- const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
- const double nJSq = calcFrobNormSq3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ);
- const double nInvJSq = calcFrobNormSqInv3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ);
- invCond(i) = 3./sqrt(nJSq*nInvJSq);
- }
- break;
- }
- }
-}
-
-void JacobianBasis::getCondNumGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &invCond) const
-{
- getCondNumGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, invCond);
-}
-
-// Calculate the inverse condition number in Frobenius norm and its gradients w.r.t. node position,
-// 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::getCondNumAndGradientsGeneral(int nJacNodes,
- const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY,
- const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals,
- fullMatrix<double> &IDI) const
-{
- fullMatrix<double> JDJ(nJacNodes, 3*numMapNodes+1);
-
- switch (_dim) {
-
- case 0 : {
- for (int i = 0; i < nJacNodes; i++) {
- for (int j = 0; j < numMapNodes; j++) {
- IDI(i,j) = 0.;
- IDI(i,j+1*numMapNodes) = 0.;
- IDI(i,j+2*numMapNodes) = 0.;
- }
- IDI(i,3*numMapNodes) = 1.;
- }
- break;
- }
-
- case 1 : {
- fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
- gSMatX.mult(nodesXYZ, dxyzdX);
- for (int i = 0; i < nJacNodes; i++) {
- const double &dxdX = dxyzdX(i,0), &dydX = dxyzdX(i,1), &dzdX = dxyzdX(i,2);
- const double &dxdY = normals(0,0), &dydY = normals(0,1), &dzdY = normals(0,2);
- const double &dxdZ = normals(1,0), &dydZ = normals(1,1), &dzdZ = normals(1,2);
- calcJDJ1D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ,
- i, numMapNodes,
- gSMatX, JDJ);
- calcIDI1D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ,
- i, numMapNodes,
- gSMatX, JDJ, IDI);
- }
- break;
- }
-
- case 2 : {
- 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);
- const double &dxdZ = normals(0,0), &dydZ = normals(0,1), &dzdZ = normals(0,2);
- calcJDJ2D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ,
- i, numMapNodes,
- gSMatX, gSMatY, JDJ);
- calcIDI2D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ,
- i, numMapNodes,
- gSMatX, gSMatY, JDJ, IDI);
- }
- break;
- }
-
- case 3 : {
- fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3), dxyzdZ(nJacNodes,3);
- 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);
- const double &dxdZ = dxyzdZ(i,0), &dydZ = dxyzdZ(i,1), &dzdZ = dxyzdZ(i,2);
- calcJDJ3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ,
- i, numMapNodes,
- gSMatX, gSMatY, gSMatZ, JDJ);
- calcIDI3D(dxdX, dxdY, dxdZ,
- dydX, dydY, dydZ,
- dzdX, dzdY, dzdZ,
- i, numMapNodes,
- gSMatX, gSMatY, gSMatZ, JDJ, IDI);
-// std::cout << "DBGTT: Jac. node " << i << ":\n";
-// std::cout << "DBGTT: -> g = " << IDI(i, 3*getNumMapNodes()) << ":\n";
-// for (int j=0; j<getNumMapNodes(); j++)
-// std::cout << "DBGTT: -> djd{x,y,z}" << j << " = (" << JDJ(i, j)
-// << ", " << JDJ(i, j+getNumMapNodes()) << ", " << JDJ(i, j+2*getNumMapNodes()) << ")\n";
-// for (int j=0; j<getNumMapNodes(); j++)
-// std::cout << "DBGTT: -> dgd{x,y,z}" << j << " = (" << IDI(i, j)
-// << ", " << IDI(i, j+getNumMapNodes()) << ", " << IDI(i, j+2*getNumMapNodes()) << ")\n";
- }
- break;
- }
- }
-}
-
-void JacobianBasis::getCondNumAndGradientsGeneral(int nJacNodes,
- const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY,
- const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals,
- fullMatrix<double> &IDI) const
-{
- getCondNumAndGradientsGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, normals, IDI);
-}
-
void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes,
const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY,
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index ace52e9907cd7474cebe56135987020b2846bea2..978e0cd1e43dd739f080e23f183441d295d95ac1 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -13,6 +13,7 @@
class GradientBasis {
public:
fullMatrix<double> gradShapeMatX, gradShapeMatY, gradShapeMatZ;
+ fullMatrix<double> gradShapeIdealMatX, gradShapeIdealMatY, gradShapeIdealMatZ;
private:
const FuncSpaceData _data;
@@ -27,6 +28,10 @@ public:
fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
fullMatrix<double> *dxyzdZ) const;
+ void getIdealGradientsFromNodes(const fullMatrix<double> &nodes,
+ fullMatrix<double> *dxyzdX,
+ fullMatrix<double> *dxyzdY,
+ fullMatrix<double> *dxyzdZ) const;
void mapFromIdealElement(fullMatrix<double> *dxyzdX,
fullMatrix<double> *dxyzdY,
fullMatrix<double> *dxyzdZ) const {
@@ -91,16 +96,6 @@ public:
getSignedIdealJacAndGradientsGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
gradShapeMatZFast,nodesXYZ,normals,JDJ);
}
- inline void getCondNumAndGradients(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullMatrix<double> &IDI) const {
- getCondNumAndGradientsGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ, nodesXYZ, normals, IDI);
- }
- inline void getCondNumAndGradientsFast(const fullMatrix<double> &nodesXYZ,
- const fullMatrix<double> &normals, fullMatrix<double> &IDI) const {
- getCondNumAndGradientsGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
- gradShapeMatZFast,nodesXYZ,normals,IDI);
- }
void getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &nodesXYZStraight,
fullVector<double> &lambdaJ , fullMatrix<double> &gradLambdaJ) const;
@@ -140,13 +135,6 @@ public:
inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
}
- inline void getCondNum(const fullMatrix<double> &nodesXYZ, fullVector<double> &invCond) const {
- getCondNumGeneral(numJacNodes, _gradBasis->gradShapeMatX,
- _gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ, nodesXYZ, invCond);
- }
- inline void getCondNumFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &invCond) const {
- getCondNumGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,invCond);
- }
//
inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
getBezier()->matrixLag2Bez.mult(jac,bez);
@@ -212,26 +200,6 @@ public:
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const;
-
- template<bool ideal>
- void getCondNumGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &invCond) const;
- void getCondNumGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, fullVector<double> &invCond) const;
-
- template<bool ideal>
- void getCondNumAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
- fullMatrix<double> &IDI) const;
- void getCondNumAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
- const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
- const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
- fullMatrix<double> &IDI) const;
-
-
};
#endif
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index 0debf3111a64cd336ddb187dd8e0060f04a80bee..f36799e280577ac86dca922fb28cc2711909dfe5 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -34,6 +34,7 @@
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "BasisFactory.h"
+#include "CondNumBasis.h"
#include "OptHomIntegralBoundaryDist.h"
#include "qualityMeasures.h"
#include "MeshOptPatch.h"
@@ -614,36 +615,38 @@ void Patch::idealJacAndGradients(int iEl, std::vector<double> &iJ, std::vector<d
}
-void Patch::initCondNum()
+void Patch::initInvCondNum()
{
// Initialize _nBezEl
- if (_nBezEl.empty()) {
- _nBezEl.resize(nEl());
- for (int iEl=0; iEl<nEl(); iEl++)
- _nBezEl[iEl] = _el[iEl]->getJacobianFuncSpace()->getNumJacNodes();
+ if (_nICNEl.empty()) {
+ _nICNEl.resize(nEl());
+ for (int iEl=0; iEl<nEl(); iEl++) {
+ const CondNumBasis *cnBasis = BasisFactory::getCondNumBasis(_el[iEl]->getTypeForMSH());
+ _nICNEl[iEl] = cnBasis->getNumCondNumNodes();
+ }
}
- // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
- // Jacobians of 3D elements
- if ((_dim == 2) && _condNormEl.empty()) {
- _condNormEl.resize(nEl());
-// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
- for (int iEl = 0; iEl < nEl(); iEl++)
- calcNormalEl2D(iEl, NS_SQRTNORM, _condNormEl[iEl], true);
- }
+// // Set normals to 2D elements (with magnitude of inverse Jacobian) or initial
+// // Jacobians of 3D elements
+// if ((_dim == 2) && _condNormEl.empty()) {
+// _condNormEl.resize(nEl());
+//// for (int iEl = 0; iEl < nEl(); iEl++) calcScaledNormalEl2D(element2entity,iEl);
+// for (int iEl = 0; iEl < nEl(); iEl++)
+// calcNormalEl2D(iEl, NS_SQRTNORM, _condNormEl[iEl], true);
+// }
}
-void Patch::condNumAndGradients(int iEl, std::vector<double> &condNum,
+void Patch::invCondNumAndGradients(int iEl, std::vector<double> &condNum,
std::vector<double> &gCondNum)
{
- const JacobianBasis *jacBasis = _el[iEl]->getJacobianFuncSpace();
- const int &numJacNodes = _nBezEl[iEl];
+ const CondNumBasis *cnBasis = BasisFactory::getCondNumBasis(_el[iEl]->getTypeForMSH());
+ const int &numICN = _nICNEl[iEl];
const int &numMapNodes = _nNodEl[iEl];
- fullMatrix<double> IDI(numJacNodes,3*numMapNodes+1);
+ fullMatrix<double> IDI(numICN, 3*numMapNodes+1);
// Coordinates of nodes
- fullMatrix<double> nodesXYZ(numMapNodes,3), normals(_dim,3);
+ fullMatrix<double> nodesXYZ(numMapNodes,3), normals;
for (int i = 0; i < numMapNodes; i++) {
int &iVi = _el2V[iEl][i];
nodesXYZ(i,0) = _xyz[iVi].x();
@@ -651,28 +654,32 @@ void Patch::condNumAndGradients(int iEl, std::vector<double> &condNum,
nodesXYZ(i,2) = _xyz[iVi].z();
}
- // Calculate Jacobian and gradients, scale if 3D (already scaled by
- // regularization normals in 2D)
- jacBasis->getCondNumAndGradients(nodesXYZ, _condNormEl[iEl], IDI);
+ // Calculate ICN and gradients
+ // TODO: Use signed measure for 2D as well
+ if (_dim == 3) {
+ cnBasis->getSignedInvCondNumAndGradients(nodesXYZ, normals, IDI);
+ }
+ else
+ cnBasis->getInvCondNumAndGradients(nodesXYZ, IDI);
// Inverse condition number
- for (int l = 0; l < numJacNodes; l++) condNum[l] = IDI(l, 3*numMapNodes);
+ for (int l = 0; l < numICN; l++) condNum[l] = IDI(l, 3*numMapNodes);
// Gradients of the inverse condition number
int iPC = 0;
- std::vector<SPoint3> gXyzV(numJacNodes);
- std::vector<SPoint3> gUvwV(numJacNodes);
+ std::vector<SPoint3> gXyzV(numICN);
+ std::vector<SPoint3> gUvwV(numICN);
for (int i = 0; i < numMapNodes; i++) {
int &iFVi = _el2FV[iEl][i];
if (iFVi >= 0) {
- for (int l = 0; l < numJacNodes; l++)
+ for (int l = 0; l < numICN; l++)
gXyzV[l] = SPoint3(IDI(l, i), IDI(l, i+numMapNodes),
IDI(l, i+2*numMapNodes));
- _coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi],gXyzV,gUvwV);
- for (int l = 0; l < numJacNodes; l++) {
- gCondNum[indGSJ(iEl,l,iPC)] = gUvwV[l][0];
- if (_nPCFV[iFVi] >= 2) gCondNum[indGSJ(iEl,l,iPC+1)] = gUvwV[l][1];
- if (_nPCFV[iFVi] == 3) gCondNum[indGSJ(iEl,l,iPC+2)] = gUvwV[l][2];
+ _coordFV[iFVi]->gXyz2gUvw(_uvw[iFVi], gXyzV, gUvwV);
+ for (int l = 0; l < numICN; l++) {
+ gCondNum[indGICN(iEl, l, iPC)] = gUvwV[l][0];
+ if (_nPCFV[iFVi] >= 2) gCondNum[indGICN(iEl, l, iPC+1)] = gUvwV[l][1];
+ if (_nPCFV[iFVi] == 3) gCondNum[indGICN(iEl, l, iPC+2)] = gUvwV[l][2];
}
iPC += _nPCFV[iFVi];
}
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index d9b477f27770d51e26f370937b0cfb4ba3c9ef7b..cc5d2bd327b27d2f75c22c77b3f7341826e9132b 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/contrib/MeshOptimizer/MeshOptPatch.h
@@ -95,11 +95,13 @@ public:
// Mesh quality
inline const int &nIJacEl(int iEl) { return _nIJacEl[iEl]; }
+ inline const int &nICNEl(int iEl) { return _nICNEl[iEl]; }
inline int indGIJac(int iEl, int l, int iPC) { return iPC*_nIJacEl[iEl]+l; }
+ inline int indGICN(int iEl, int l, int iPC) { return iPC*_nICNEl[iEl]+l; }
void initIdealJac();
void idealJacAndGradients(int iEl, std::vector<double> &NCJ, std::vector<double> &gNCJ);
- void initCondNum();
- void condNumAndGradients(int iEl, std::vector<double> &condNum, std::vector<double> &gCondNum);
+ void initInvCondNum();
+ void invCondNumAndGradients(int iEl, std::vector<double> &condNum, std::vector<double> &gCondNum);
private:
@@ -155,7 +157,8 @@ private:
std::vector<int> _nIJacEl; // Number of NCJ values for an el
std::vector<fullMatrix<double> > _IJacNormEl; // Normals to 2D elements for Jacobian regularization and scaling
std::vector<double> _invIJac; // Initial Jacobians for 3D elements
- std::vector<fullMatrix<double> > _condNormEl; // Normals to 2D elements for inverse conditioning computation
+ std::vector<int> _nICNEl; // Number of inv. cond. number values for an el.
+// std::vector<fullMatrix<double> > _condNormEl; // Normals to 2D elements for inverse conditioning computation
};
diff --git a/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h b/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
index 65b758fb3f0d91c9f8e5daa85bea529845d39a6d..3ef5bde3c2bcf7fb5bf87dfbc0e694c9429c0636 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
+++ b/contrib/MeshQualityOptimizer/MeshQualityObjContribInvCond.h
@@ -8,11 +8,11 @@
template<class FuncType>
-class ObjContribInvCond : public ObjContrib, public FuncType
+class ObjContribInvCondNum : public ObjContrib, public FuncType
{
public:
- ObjContribInvCond(double weight);
- virtual ~ObjContribInvCond() {}
+ ObjContribInvCondNum(double weight);
+ virtual ~ObjContribInvCondNum() {}
virtual ObjContrib *copy() const;
virtual void initialize(Patch *mesh);
virtual bool fail() { return false; }
@@ -29,45 +29,45 @@ protected:
template<class FuncType>
-ObjContribInvCond<FuncType>::ObjContribInvCond(double weight) :
- ObjContrib("InvCond", FuncType::getNamePrefix()+"InvCond"),
+ObjContribInvCondNum<FuncType>::ObjContribInvCondNum(double weight) :
+ ObjContrib("InvCondNum", FuncType::getNamePrefix()+"InvCondNum"),
_mesh(0), _weight(weight)
{
}
template<class FuncType>
-ObjContrib *ObjContribInvCond<FuncType>::copy() const
+ObjContrib *ObjContribInvCondNum<FuncType>::copy() const
{
- return new ObjContribInvCond<FuncType>(*this);
+ return new ObjContribInvCondNum<FuncType>(*this);
}
template<class FuncType>
-void ObjContribInvCond<FuncType>::initialize(Patch *mesh)
+void ObjContribInvCondNum<FuncType>::initialize(Patch *mesh)
{
_mesh = mesh;
- _mesh->initCondNum();
+ _mesh->initInvCondNum();
updateMinMax();
FuncType::initialize(_min, _max);
}
template<class FuncType>
-bool ObjContribInvCond<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
+bool ObjContribInvCondNum<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
{
_min = BIGVAL;
_max = -BIGVAL;
for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
- std::vector<double> invCond(_mesh->nBezEl(iEl)); // Min. of Metric
- std::vector<double> gInvCond(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
- _mesh->condNumAndGradients(iEl, invCond, gInvCond);
- for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Add contribution for each Bezier coeff.
+ std::vector<double> invCond(_mesh->nICNEl(iEl)); // Min. of Metric
+ std::vector<double> gInvCond(_mesh->nICNEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
+ _mesh->invCondNumAndGradients(iEl, invCond, gInvCond);
+ for (int l = 0; l < _mesh->nICNEl(iEl); l++) { // Add contribution for each Bezier coeff.
Obj += _weight * FuncType::compute(invCond[l]);
const double dfact = _weight * FuncType::computeDiff(invCond[l]);
for (int iPC = 0; iPC < _mesh->nPCEl(iEl); iPC++)
- gradObj[_mesh->indPCEl(iEl,iPC)] += dfact * gInvCond[_mesh->indGSJ(iEl, l, iPC)];
+ gradObj[_mesh->indPCEl(iEl,iPC)] += dfact * gInvCond[_mesh->indGICN(iEl, l, iPC)];
_min = std::min(_min, invCond[l]);
_max = std::max(_max, invCond[l]);
}
@@ -78,16 +78,16 @@ bool ObjContribInvCond<FuncType>::addContrib(double &Obj, alglib::real_1d_array
template<class FuncType>
-void ObjContribInvCond<FuncType>::updateMinMax()
+void ObjContribInvCondNum<FuncType>::updateMinMax()
{
_min = BIGVAL;
_max = -BIGVAL;
for (int iEl = 0; iEl < _mesh->nEl(); iEl++) {
- std::vector<double> invCond(_mesh->nBezEl(iEl)); // Min. of Metric
- std::vector<double> dumGInvCond(_mesh->nBezEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
- _mesh->condNumAndGradients(iEl, invCond, dumGInvCond);
- for (int l = 0; l < _mesh->nBezEl(iEl); l++) { // Add contribution for each Bezier coeff.
+ std::vector<double> invCond(_mesh->nICNEl(iEl)); // Min. of Metric
+ std::vector<double> dumGInvCond(_mesh->nICNEl(iEl)*_mesh->nPCEl(iEl)); // Dummy gradients of metric min.
+ _mesh->invCondNumAndGradients(iEl, invCond, dumGInvCond);
+ for (int l = 0; l < _mesh->nICNEl(iEl); l++) { // Add contribution for each Bezier coeff.
_min = std::min(_min, invCond[l]);
_max = std::max(_max, invCond[l]);
}
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
index db96312c383c825b9d9800d987aaa2737a6c4934..8c47078ea457c9e38562c2c7fe11d029f5f483d1 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.cpp
@@ -26,7 +26,7 @@ struct QualPatchDefParameters : public MeshOptPatchDef
double limDist, MElement *el) const;
private:
bool _excludeHex, _excludePrism;
- double _idealJacMin;
+ double _idealJacMin, _invCondNumMin;
double _distanceFactor;
};
@@ -36,6 +36,7 @@ QualPatchDefParameters::QualPatchDefParameters(const MeshQualOptParameters &p)
_excludeHex = p.excludeHex;
_excludePrism = p.excludePrism;
_idealJacMin = p.minTargetIdealJac;
+ _invCondNumMin = p.minTargetInvCondNum;
strategy = (p.strategy == 1) ? MeshOptParameters::STRAT_ONEBYONE :
MeshOptParameters::STRAT_CONNECTED;
minLayers = (p.dim == 3) ? 1 : 0;
@@ -56,9 +57,12 @@ double QualPatchDefParameters::elBadness(MElement *el) const
const int typ = el->getType();
if (_excludeHex && (typ == TYPE_HEX)) return 1.;
if (_excludePrism && (typ == TYPE_PRI)) return 1.;
- double jMin, jMax;
- el->idealJacRange(jMin, jMax);
- return jMin-_idealJacMin;
+// double jMin, jMax;
+// el->idealJacRange(jMin, jMax);
+// return jMin-_idealJacMin;
+ double iCNMin, iCNMax;
+ el->invCondNumRange(iCNMin, iCNMax);
+ return iCNMin-_invCondNumMin;
}
@@ -94,6 +98,8 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
Patch::LS_MINEDGELENGTH);
ObjContribIdealJac<ObjContribFuncBarrierMovMin> minIdealJacBarFunc(1.);
minIdealJacBarFunc.setTarget(p.minTargetIdealJac, 1.);
+ ObjContribInvCondNum<ObjContribFuncBarrierMovMin> minInvCondNumBarFunc(1.);
+ minInvCondNumBarFunc.setTarget(p.minTargetInvCondNum, 1.);
MeshOptPass minJacPass;
minJacPass.barrierIterMax = p.optPassMax;
@@ -102,9 +108,18 @@ void MeshQualityOptimizer(GModel *gm, MeshQualOptParameters &p)
minJacPass.contrib.push_back(&minIdealJacBarFunc);
par.pass.push_back(minJacPass);
+ MeshOptPass minInvCondNumPass;
+ minInvCondNumPass.barrierIterMax = p.optPassMax;
+ minInvCondNumPass.optIterMax = p.itMax;
+ minInvCondNumPass.contrib.push_back(&nodeDistFunc);
+ minInvCondNumPass.contrib.push_back(&minInvCondNumBarFunc);
+ par.pass.push_back(minInvCondNumPass);
+
meshOptimizer(gm, par);
p.CPU = par.CPU;
- p.minIdealJac = minIdealJacBarFunc.getMin();
- p.maxIdealJac = minIdealJacBarFunc.getMax();
+// p.minIdealJac = minIdealJacBarFunc.getMin();
+// p.maxIdealJac = minIdealJacBarFunc.getMax();
+ p.minInvCondNum = minInvCondNumBarFunc.getMin();
+ p.maxInvCondNum = minInvCondNumBarFunc.getMax();
}
diff --git a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
index 36722f0a453d634a96f480992703533d7841cbe8..052addc17207b99cd4b77e0d30e0be39453c36c9 100644
--- a/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
+++ b/contrib/MeshQualityOptimizer/MeshQualityOptimizer.h
@@ -35,6 +35,7 @@ class GModel;
struct MeshQualOptParameters {
bool excludeHex, excludePrism;
double minTargetIdealJac;
+ double minTargetInvCondNum;
double weightFixed; // weight of the energy for fixed nodes
double weightFree; // weight of the energy for free nodes
int nbLayers; // number of layers taken around a bad element
@@ -52,15 +53,18 @@ struct MeshQualOptParameters {
double adaptBlobDistFact; // Growth factor in distance factor for blob adaptation
int SUCCESS ; // 0 --> success , 1 --> Not converged
- double minIdealJac, maxIdealJac; // after optimization, range of jacobians
+// double minIdealJac, maxIdealJac; // after optimization, range of jacobians
+ double minInvCondNum, maxInvCondNum; // after optimization, range of jacobians
double CPU; // Time for optimization
MeshQualOptParameters ()
- : excludeHex(false), excludePrism(false), minTargetIdealJac(0.1), weightFixed(1000.),
+ : excludeHex(false), excludePrism(false), minTargetIdealJac(0.1),
+ minTargetInvCondNum(0.1), weightFixed(1000.),
weightFree (1.), nbLayers (6) , dim(3) , itMax(300), optPassMax(50),
onlyVisible(true), distanceFactor(12), fixBndNodes(false), strategy(0),
maxAdaptBlob(3), adaptBlobLayerFact(2.), adaptBlobDistFact(2.), CPU(0.),
- minIdealJac(0.), maxIdealJac(0.), SUCCESS(-1)
+// minIdealJac(0.), maxIdealJac(0.), SUCCESS(-1)
+ minInvCondNum(0.), maxInvCondNum(0.), SUCCESS(-1)
{
}
};