From a3c2a6c92f6cb9e7be4ee775a6e10c30b6751eaf Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Wed, 17 Sep 2014 14:02:20 +0000
Subject: [PATCH] Disabled mapping from ideal element when calculatin Jacobian,
introduced in r19425
---
Numeric/JacobianBasis.cpp | 50 +++++++++++++++++++++++++--------------
Numeric/JacobianBasis.h | 9 +++++--
Numeric/MetricBasis.cpp | 11 +++++----
Numeric/MetricBasis.h | 1 +
4 files changed, 46 insertions(+), 25 deletions(-)
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 7b8980fdd9..485910cae1 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -322,7 +322,7 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
// 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>
+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
@@ -362,7 +362,7 @@ 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 (!scaling) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, NULL);
+ 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);
@@ -378,7 +378,7 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
gSMatZ.mult(nodesXYZ, dxyzdZ);
- if (!scaling) _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, &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);
@@ -402,7 +402,7 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
- getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
+ getJacobianGeneral<false, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
}
// Calculate (signed) scaled Jacobian for one element, with normal vectors to straight element
@@ -411,13 +411,13 @@ void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<dou
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
- getJacobianGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
+ getJacobianGeneral<true, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, 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>
+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,
@@ -465,7 +465,7 @@ 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 (!scaling) {
+ if (ideal) {
_gradBasis->mapFromIdealElement(&dxdX, &dxdY, NULL);
_gradBasis->mapFromIdealElement(&dydX, &dydY, NULL);
_gradBasis->mapFromIdealElement(&dzdX, &dzdY, NULL);
@@ -500,7 +500,7 @@ 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 (!scaling) {
+ if (ideal) {
_gradBasis->mapFromIdealElement(&dxdX, &dxdY, &dxdZ);
_gradBasis->mapFromIdealElement(&dydX, &dydY, &dydZ);
_gradBasis->mapFromIdealElement(&dzdX, &dzdY, &dzdZ);
@@ -535,7 +535,7 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
- getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
+ getJacobianGeneral<false, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
}
// Calculate (signed) scaled Jacobian for several elements, with normal vectors to straight elements
@@ -545,17 +545,19 @@ void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<dou
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
- getJacobianGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
+ getJacobianGeneral<true, false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, 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.
-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
+template<bool ideal>
+inline 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
{
switch (_dim) {
@@ -593,7 +595,7 @@ void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes, const fullMat
fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
- _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, NULL);
+ 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);
@@ -623,7 +625,7 @@ void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes, const fullMat
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
gSMatZ.mult(nodesXYZ, dxyzdZ);
- _gradBasis->mapFromIdealElement(&dxyzdX, &dxyzdY, &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);
@@ -654,6 +656,18 @@ void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes, const fullMat
}
+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::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
const fullMatrix<double> &nodesXYZStraight,
fullVector<double> &lambdaJ, fullMatrix<double> &gradLambdaJ) const
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index ca744f6e57..4dd280ecc8 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -125,11 +125,11 @@ class JacobianBasis {
private :
- template<bool scaling>
+ template<bool scaling, bool ideal>
inline 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>
+ template<bool scaling, bool ideal>
inline void getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
@@ -149,6 +149,11 @@ class JacobianBasis {
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
+ template<bool ideal>
+ inline 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,
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index 189a5d2953..36cbfbb9db 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -95,7 +95,7 @@ double MetricBasis::minRCorner(MElement *el)
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff(el->getDim(), gradients, nodes, metCoeffLag);
+ _fillCoeff<false>(el->getDim(), gradients, nodes, metCoeffLag);
// Compute min_corner(R)
return _computeMinlagR(jacLag, metCoeffLag, nSampPnts);
@@ -183,7 +183,7 @@ double MetricBasis::getBoundMinR(MElement *el)
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff(el->getDim(), _gradients, nodes, metCoeffLag);
+ _fillCoeff<false>(el->getDim(), _gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
@@ -808,7 +808,7 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
// Metric coefficients
fullMatrix<double> metCoeffLag;
- _fillCoeff(el->getDim(), _gradients, nodes, metCoeffLag);
+ _fillCoeff<false>(el->getDim(), _gradients, nodes, metCoeffLag);
fullMatrix<double> *metCoeff;
metCoeff = new fullMatrix<double>(nSampPnts, metCoeffLag.size2());
_bezier->matrixLag2Bez.mult(metCoeffLag, *metCoeff);
@@ -816,6 +816,7 @@ void MetricBasis::_getMetricData(MElement *el, MetricData *&md) const
md = new MetricData(metCoeff, jac, -1, 0, 0);
}
+template<bool ideal>
void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
fullMatrix<double> &nodes, fullMatrix<double> &coeff)
{
@@ -832,7 +833,7 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
{
fullMatrix<double> dxydX(nSampPnts,3), dxydY(nSampPnts,3);
gradients->getGradientsFromNodes(nodes, &dxydX, &dxydY, NULL);
- gradients->mapFromIdealElement(&dxydX, &dxydY, NULL);
+ if (ideal) gradients->mapFromIdealElement(&dxydX, &dxydY, NULL);
coeff.resize(nSampPnts, 3);
for (int i = 0; i < nSampPnts; i++) {
@@ -851,7 +852,7 @@ void MetricBasis::_fillCoeff(int dim, const GradientBasis *gradients,
{
fullMatrix<double> dxyzdX(nSampPnts,3), dxyzdY(nSampPnts,3), dxyzdZ(nSampPnts,3);
gradients->getGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
- gradients->mapFromIdealElement(&dxyzdX, &dxyzdY, &dxyzdZ);
+ if (ideal) gradients->mapFromIdealElement(&dxyzdX, &dxyzdY, &dxyzdZ);
coeff.resize(nSampPnts, 7);
for (int i = 0; i < nSampPnts; i++) {
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index bab7253e53..fc15f45a5c 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -78,6 +78,7 @@ private:
void _getMetricData(MElement*, MetricData*&) const;
double _subdivideForRmin(MetricData*, double RminLag, double tol) const;
+ template<bool ideal>
static void _fillCoeff(int dim, const GradientBasis*,
fullMatrix<double> &nodes, fullMatrix<double> &coeff);
static double _computeMinlagR(const fullVector<double> &jac,
--
GitLab