From 66fa13429593ddcc28f1a92c8706f51812cc963f Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Fri, 17 Oct 2014 14:49:07 +0000
Subject: [PATCH] add faster and more accurate function
MetricBasis::minSampledRnew(..)
---
Numeric/BasisFactory.h | 3 +
Numeric/FuncSpaceData.cpp | 9 +--
Numeric/FuncSpaceData.h | 10 ++--
Numeric/JacobianBasis.cpp | 32 ++++++++++
Numeric/JacobianBasis.h | 1 +
Numeric/MetricBasis.cpp | 123 ++++++++++++++++++++++++++++----------
Numeric/MetricBasis.h | 1 +
Numeric/bezierBasis.cpp | 21 +++++--
Numeric/bezierBasis.h | 3 +
9 files changed, 161 insertions(+), 42 deletions(-)
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index d5c509c164..b1a11a4c34 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -61,6 +61,9 @@ class BasisFactory
static const GradientBasis* getGradientBasis(int tag, int order) {
return getGradientBasis(FuncSpaceData(true, tag, order));
}
+ static const GradientBasis* getGradientBasis(int tag) {
+ return getGradientBasis(FuncSpaceData(tag));
+ }
// Bezier
static const bezierBasis* getBezierBasis(FuncSpaceData);
diff --git a/Numeric/FuncSpaceData.cpp b/Numeric/FuncSpaceData.cpp
index 8e0dfc9382..586213617e 100644
--- a/Numeric/FuncSpaceData.cpp
+++ b/Numeric/FuncSpaceData.cpp
@@ -6,19 +6,20 @@
#include "FuncSpaceData.h"
#include "MElement.h"
-FuncSpaceData::FuncSpaceData(MElement *el, int *serendip) :
+FuncSpaceData::FuncSpaceData(const MElement *el, const bool *serendip) :
_tag(el->getTypeForMSH()), _spaceOrder(el->getPolynomialOrder()),
_nij(0), _nk(_spaceOrder), _pyramidalSpace(el->getType() == TYPE_PYR),
_serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
{}
-FuncSpaceData::FuncSpaceData(MElement *el, int order, int *serendip) :
+FuncSpaceData::FuncSpaceData(const MElement *el, int order, const bool *serendip) :
_tag(el->getTypeForMSH()), _spaceOrder(order),
_nij(0), _nk(_spaceOrder), _pyramidalSpace(el->getType() == TYPE_PYR),
_serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
{}
-FuncSpaceData::FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *serendip) :
+FuncSpaceData::FuncSpaceData(const MElement *el, bool pyr, int nij, int nk,
+ const bool *serendip) :
_tag(el->getTypeForMSH()), _spaceOrder(pyr ? nij+nk : std::max(nij, nk)),
_nij(nij), _nk(nk), _pyramidalSpace(pyr),
_serendipity(serendip ? *serendip : el->getIsOnlySerendipity())
@@ -27,7 +28,7 @@ FuncSpaceData::FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *seren
Msg::Error("Creation of pyramidal space data for a non-pyramid element !");
}
-FuncSpaceData::FuncSpaceData(int tag, int *serendip) :
+FuncSpaceData::FuncSpaceData(int tag, const bool *serendip) :
_tag(tag), _spaceOrder(ElementType::OrderFromTag(tag)),
_nij(0), _nk(_spaceOrder),
_pyramidalSpace(ElementType::ParentTypeFromTag(tag) == TYPE_PYR),
diff --git a/Numeric/FuncSpaceData.h b/Numeric/FuncSpaceData.h
index 0ddc82b565..649a5b9413 100644
--- a/Numeric/FuncSpaceData.h
+++ b/Numeric/FuncSpaceData.h
@@ -31,12 +31,14 @@ public:
_pyramidalSpace(false), _serendipity(false) {}
// Constructors using MElement*
- FuncSpaceData(MElement *el, int *serendip = NULL);
- FuncSpaceData(MElement *el, int order, int *serendip = NULL);
- FuncSpaceData(MElement *el, bool pyr, int nij, int nk, int *serendip = NULL);
+ FuncSpaceData(const MElement *el, const bool *serendip = NULL);
+ FuncSpaceData(const MElement *el, int order, const bool *serendip = NULL);
+ FuncSpaceData(const MElement *el,
+ bool pyr, int nij, int nk,
+ const bool *serendip = NULL);
// Constructors using element tag
- FuncSpaceData(int tag, int *serendip = NULL);
+ FuncSpaceData(int tag, const bool *serendip = NULL);
// constructors using element tag or element type
FuncSpaceData(bool isTag, int tagOrType, int order,
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index 2bf7bc4fe0..a874038d32 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -854,3 +854,35 @@ int JacobianBasis::jacobianOrder(int parentType, int order)
return 0;
}
}
+
+FuncSpaceData JacobianBasis::jacobianMatrixSpace(int type, int order)
+{
+ if (type == TYPE_PYR) {
+ Msg::Fatal("jacobianMatrixSpace not yet implemented for pyramids");
+ return FuncSpaceData(false, type, false, 1, 0);
+ }
+ int jacOrder = -1;
+ switch (type) {
+ case TYPE_PNT :
+ jacOrder = 0;
+ break;
+
+ case TYPE_LIN :
+ case TYPE_TRI :
+ case TYPE_TET :
+ jacOrder = order - 1;
+ break;
+
+ case TYPE_QUA :
+ case TYPE_PRI :
+ case TYPE_HEX :
+ jacOrder = order;
+ break;
+
+ default :
+ Msg::Error("Unknown element type %d, return order 0", type);
+ return 0;
+ }
+
+ return FuncSpaceData(true, ElementType::getTag(type, order), jacOrder);
+}
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 7bfbf85375..bbcca1b059 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -202,6 +202,7 @@ public:
//
static int jacobianOrder(int tag);
static int jacobianOrder(int parentType, int order);
+ static FuncSpaceData jacobianMatrixSpace(int type, int order);
private :
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index eb4959eec4..877e9c4c1e 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -153,36 +153,8 @@ double MetricBasis::minSampledR(MElement *el, int order)
double MetricBasis::getMinSampledR(MElement *el, int deg) const
{
fullMatrix<double> samplingPoints;
-
- switch (el->getType()) {
- case TYPE_PNT :
- samplingPoints = gmshGeneratePointsLine(0);
- break;
- case TYPE_LIN :
- samplingPoints = gmshGeneratePointsLine(deg);
- break;
- case TYPE_TRI :
- samplingPoints = gmshGeneratePointsTriangle(deg,false);
- break;
- case TYPE_QUA :
- samplingPoints = gmshGeneratePointsQuadrangle(deg,false);
- break;
- case TYPE_TET :
- samplingPoints = gmshGeneratePointsTetrahedron(deg,false);
- break;
- case TYPE_PRI :
- samplingPoints = gmshGeneratePointsPrism(deg,false);
- break;
- case TYPE_HEX :
- samplingPoints = gmshGeneratePointsHexahedron(deg,false);
- break;
- case TYPE_PYR :
- samplingPoints = gmshGeneratePointsPyramidGeneral(true, deg+1, deg);
- break;
- default :
- Msg::Error("Unknown Jacobian function space for element type %d", el->getType());
- return -1;
- }
+ bool serendip = false;
+ gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), samplingPoints);
MetricData *md;
_getMetricData(el, md);
@@ -199,6 +171,73 @@ double MetricBasis::getMinSampledR(MElement *el, int deg) const
return min;
}
+double MetricBasis::minSampledRnew(MElement *el, int deg)
+{
+ int dim = el->getDim();
+ if (dim < 2) return 1;
+
+ fullMatrix<double> nodes(el->getNumVertices(), 3);
+ el->getNodesCoord(nodes);
+
+ // if deg is small: sample at every point
+ // otherwise: compute Bezier coefficients and interpolate
+ //const bool BezierInterp = deg > 100;
+
+ FuncSpaceData jacMatSpace;
+ /*if (BezierInterp) {
+ jacMatSpace =
+ JacobianBasis::jacobianMatrixSpace(el->getType(), el->getPolynomialOrder());
+ }
+ else {*/
+ const bool serendip = false;
+ jacMatSpace = FuncSpaceData(el, deg, &serendip);
+ //}
+
+ const GradientBasis *gb = BasisFactory::getGradientBasis(jacMatSpace);
+ fullMatrix<double> coeffJacMat(gb->getNumSamplingPoints(), 3*dim);
+ fullMatrix<double> dxyzdX(coeffJacMat, 0, 3);
+ fullMatrix<double> dxyzdY(coeffJacMat, 3, 3);
+ if (dim == 2)
+ gb->getIdealGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, NULL);
+ else {
+ fullMatrix<double> dxyzdZ(coeffJacMat, 6, 3);
+ gb->getIdealGradientsFromNodes(nodes, &dxyzdX, &dxyzdY, &dxyzdZ);
+ }
+
+ /*if (BezierInterp) {
+ const bezierBasis *bb = BasisFactory::getBezierBasis(jacMatSpace);
+ fullMatrix<double> points;
+ const bool serendip = false;
+ gmshGeneratePoints(FuncSpaceData(el, deg, &serendip), points);
+
+ fullMatrix<double> bezCoeffJacMat;
+ bb->lag2Bez(coeffJacMat, bezCoeffJacMat);
+ bb->interpolate(bezCoeffJacMat, points, coeffJacMat);
+ }*/
+
+ double minR = 1;
+ for (int k = 0; k < coeffJacMat.size1(); ++k) {
+ fullMatrix<double> metric(dim, dim);
+ for (int i = 0; i < dim; ++i) {
+ for (int j = i; j < dim; ++j) {
+ for (int n = 0; n < 3; ++n)
+ metric(i, j) += coeffJacMat(k, 3*i+n) * coeffJacMat(k, 3*j+n);
+ }
+ }
+ metric(1, 0) = metric(0, 1);
+ if (dim == 3) {
+ metric(2, 0) = metric(0, 2);
+ metric(2, 1) = metric(1, 2);
+ }
+ fullVector<double> valReal(dim), valImag(dim);
+ fullMatrix<double> vecLeft(dim, dim), vecRight(dim, dim);
+ metric.eig(valReal, valImag, vecLeft, vecRight, true);
+ minR = std::min(minR, std::sqrt(valReal(0) / valReal(dim-1)));
+ }
+
+ return minR;
+}
+
double MetricBasis::getBoundMinR(MElement *el) const
{
int nSampPnts = _gradients->getNumSamplingPoints();
@@ -607,6 +646,30 @@ bool MetricBasis::validateBezierForMetricAndJacobian()
metricBasis->interpolateAfterNSubdivisions(el, numSubdiv, numSampPnt,
isub, uvw, metric_Bez);
+ /*{
+ static int deg = 2;
+ MetricBasis::minSampledR(el, deg);
+ MetricBasis::minSampledRnew(el, deg);
+
+ double time = Cpu();
+ double minR1;
+ for (int cn = 0; cn < 100; ++cn)
+ minR1 = MetricBasis::minSampledR(el, deg);
+ double tm1 = Cpu()-time;
+
+ time = Cpu();
+ double minR2;
+ for (int cn = 0; cn < 100; ++cn)
+ minR2 = MetricBasis::minSampledRnew(el, deg);
+ double tm2 = Cpu()-time;
+
+ if (std::abs(minR1-minR2) < 1e-14)
+ Msg::Info("ok deg %d, times %g %g speedup %g", deg, tm1, tm2, tm1/tm2);
+ else
+ Msg::Error("not ok deg %d: %g vs %g, times %g %g speedup %g", deg, minR1, minR2, tm1, tm2, tm1/tm2);
+ //++deg;
+ }*/
+
int numBadMatch = 0;
int numBadMatchTensor = 0;
double maxBadMatch = 0;
diff --git a/Numeric/MetricBasis.h b/Numeric/MetricBasis.h
index 2416f8d12f..0817458ace 100644
--- a/Numeric/MetricBasis.h
+++ b/Numeric/MetricBasis.h
@@ -64,6 +64,7 @@ public:
static double boundMinR(MElement *el);
static double minSampledR(MElement *el, int order);
+ static double minSampledRnew(MElement *el, int order);
double getBoundMinR(MElement*) const;
double getMinSampledR(MElement*, int order) const;
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 5eefae9b27..2a56445b35 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -499,7 +499,7 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
fullMatrix<double> &result,
bool bezCoord) const
{
- if (result.size1() < uvw.size1() || result.size2() < coeffs.size2())
+ if (result.size1() != uvw.size1() || result.size2() != coeffs.size2())
result.resize(uvw.size1(), coeffs.size2());
fullMatrix<double> bezuvw = uvw;
@@ -535,11 +535,24 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
}
}
+void bezierBasis::lag2Bez(const fullMatrix<double> &lag,
+ fullMatrix<double> &bez) const
+{
+ if (lag.size1() != matrixLag2Bez.size1()) {
+ Msg::Error("matrix not the right size in lag2Bez function %d vs %d",
+ lag.size1(), matrixLag2Bez.size1());
+ }
+ if (bez.size1() != lag.size1() || bez.size2() != lag.size2()) {
+ bez.resize(lag.size1(), lag.size2());
+ }
+ matrixLag2Bez.mult(lag, bez);
+}
+
void bezierBasis::subdivideBezCoeff(const fullMatrix<double> &coeff,
fullMatrix<double> &subCoeff) const
{
- if (subCoeff.size1() < subDivisor.size1()
- || subCoeff.size2() < coeff.size2() ) {
+ if (subCoeff.size1() != subDivisor.size1()
+ || subCoeff.size2() != coeff.size2() ) {
subCoeff.resize(subDivisor.size1(), coeff.size2());
}
subDivisor.mult(coeff, subCoeff);
@@ -548,7 +561,7 @@ void bezierBasis::subdivideBezCoeff(const fullMatrix<double> &coeff,
void bezierBasis::subdivideBezCoeff(const fullVector<double> &coeff,
fullVector<double> &subCoeff) const
{
- if (subCoeff.size() < subDivisor.size1()) {
+ if (subCoeff.size() != subDivisor.size1()) {
subCoeff.resize(subDivisor.size1());
}
subDivisor.mult(coeff, subCoeff);
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index f08e305a9e..b301391fd2 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -48,6 +48,9 @@ class bezierBasis {
// generate Bezier points
void generateBezierPoints(fullMatrix<double>&) const;
+ // transform coeff Lagrange into Bezier coeff
+ void lag2Bez(const fullMatrix<double> &lag, fullMatrix<double> &bez) const;
+
// Subdivide Bezier coefficients
void subdivideBezCoeff(const fullMatrix<double> &coeff,
fullMatrix<double> &subCoeff) const;
--
GitLab