From 542bdc6f1deea0c526b25e2ac364db7a0f57a666 Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Thu, 31 Jul 2014 12:35:28 +0000
Subject: [PATCH] (1) Faster JacobianBasis construction: construct bezierBasis
only if necessary (2) Make pyramid's points and monomials ready for
differenciation of subdivision-ready and optimization-ready basis (3)
Correction number lagragian bezier coeff for pyramids p1
---
Numeric/BasisFactory.h | 9 +-
Numeric/JacobianBasis.cpp | 62 +++++++------
Numeric/JacobianBasis.h | 19 ++--
Numeric/MetricBasis.cpp | 16 ++--
Numeric/bezierBasis.cpp | 106 ++++++++++-----------
Numeric/bezierBasis.h | 20 ++--
Numeric/pointsGenerators.cpp | 175 +++++++++++++++++++++++++++++++++--
Numeric/pointsGenerators.h | 9 +-
Plugin/AnalyseCurvedMesh.cpp | 25 ++---
Plugin/AnalyseCurvedMesh.h | 3 -
10 files changed, 305 insertions(+), 139 deletions(-)
diff --git a/Numeric/BasisFactory.h b/Numeric/BasisFactory.h
index 6c4e4e110c..82b5c3731a 100644
--- a/Numeric/BasisFactory.h
+++ b/Numeric/BasisFactory.h
@@ -9,6 +9,7 @@
#include "MElement.h"
#include "MPyramid.h"
#include "nodalBasis.h"
+#include "bezierBasis.h"
#include "JacobianBasis.h"
#include "MetricBasis.h"
@@ -33,10 +34,10 @@ class BasisFactory
static const GradientBasis* getGradientBasis(int tag, int order);
static const bezierBasis* getBezierBasis(int parentTag, int order);
- static inline const bezierBasis* getBezierBasis(int tag) {
- return getBezierBasis(ElementType::ParentTypeFromTag(tag),
- ElementType::OrderFromTag(tag) );
- }
+ static const bezierBasis* getBezierBasis(int tag) {
+ return getBezierBasis(ElementType::ParentTypeFromTag(tag),
+ ElementType::OrderFromTag(tag) );
+ }
static void clearAll();
};
diff --git a/Numeric/JacobianBasis.cpp b/Numeric/JacobianBasis.cpp
index ba641a7ccc..0d1b868377 100644
--- a/Numeric/JacobianBasis.cpp
+++ b/Numeric/JacobianBasis.cpp
@@ -3,6 +3,8 @@
// See the LICENSE.txt file for license information. Please report all
// bugs and problems to the public mailing list <gmsh@geuz.org>.
+#include "JacobianBasis.h"
+
#include "GmshDefines.h"
#include "GmshMessage.h"
@@ -88,15 +90,11 @@ void GradientBasis::getGradientsFromNodes(const fullMatrix<double> &nodes,
if (dxyzdZ) gradShapeMatZ.mult(nodes, *dxyzdZ);
}
-JacobianBasis::JacobianBasis(int tag, int jacOrder)
+JacobianBasis::JacobianBasis(int tag, int jacOrder) :
+ _bezier(NULL), _tag(tag), _dim(ElementType::DimensionFromTag(tag)),
+ _jacOrder(jacOrder >= 0 ? jacOrder : JacobianBasis::jacobianOrder(tag))
{
const int parentType = ElementType::ParentTypeFromTag(tag);
- const int order = ElementType::OrderFromTag(tag);
- int jacobianOrder;
- if (jacOrder < 0)
- jacobianOrder = JacobianBasis::jacobianOrder(parentType, order);
- else
- jacobianOrder = jacOrder;
const int primJacobianOrder = JacobianBasis::jacobianOrder(parentType, 1);
fullMatrix<double> lagPoints; // Sampling points
@@ -106,25 +104,25 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder)
lagPoints = gmshGeneratePointsLine(0);
break;
case TYPE_LIN :
- lagPoints = gmshGeneratePointsLine(jacobianOrder);
+ lagPoints = gmshGeneratePointsLine(_jacOrder);
break;
case TYPE_TRI :
- lagPoints = gmshGeneratePointsTriangle(jacobianOrder,false);
+ lagPoints = gmshGeneratePointsTriangle(_jacOrder,false);
break;
case TYPE_QUA :
- lagPoints = gmshGeneratePointsQuadrangle(jacobianOrder,false);
+ lagPoints = gmshGeneratePointsQuadrangle(_jacOrder,false);
break;
case TYPE_TET :
- lagPoints = gmshGeneratePointsTetrahedron(jacobianOrder,false);
+ lagPoints = gmshGeneratePointsTetrahedron(_jacOrder,false);
break;
case TYPE_PRI :
- lagPoints = gmshGeneratePointsPrism(jacobianOrder,false);
+ lagPoints = gmshGeneratePointsPrism(_jacOrder,false);
break;
case TYPE_HEX :
- lagPoints = gmshGeneratePointsHexahedron(jacobianOrder,false);
+ lagPoints = gmshGeneratePointsHexahedron(_jacOrder,false);
break;
case TYPE_PYR :
- lagPoints = generateJacPointsPyramid(jacobianOrder);
+ lagPoints = generateJacPointsPyramid(_jacOrder);
break;
default :
Msg::Error("Unknown Jacobian function space for element tag %d", tag);
@@ -132,11 +130,8 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder)
}
numJacNodes = lagPoints.size1();
- // Store Bezier basis
- bezier = BasisFactory::getBezierBasis(parentType, jacobianOrder);
-
// Store shape function gradients of mapping at Jacobian nodes
- _gradBasis = BasisFactory::getGradientBasis(tag, jacobianOrder);
+ _gradBasis = BasisFactory::getGradientBasis(tag, _jacOrder);
// Compute matrix for lifting from primary Jacobian basis to Jacobian basis
int primJacType = ElementType::getTag(parentType, primJacobianOrder, false);
@@ -203,6 +198,14 @@ JacobianBasis::JacobianBasis(int tag, int jacOrder)
}
}
+const bezierBasis* JacobianBasis::getBezier() const {
+ if (_bezier) return _bezier;
+ const int parentType = ElementType::ParentTypeFromTag(_tag);
+ const_cast<JacobianBasis*>(this)->_bezier =
+ BasisFactory::getBezierBasis(parentType, _jacOrder);
+ return _bezier;
+}
+
// Computes (unit) normals to straight line element at barycenter (with norm of gradient as return value)
double JacobianBasis::getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const
{
@@ -286,10 +289,7 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
-
- const int dim = bezier->getDim();
-
- switch (dim) {
+ switch (_dim) {
case 0 : {
for (int i = 0; i < nJacNodes; i++) jacobian(i) = 1.;
@@ -383,8 +383,7 @@ inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<do
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
-
- switch (bezier->getDim()) {
+ switch (_dim) {
case 0 : {
const int numEl = nodesX.size2();
@@ -508,8 +507,7 @@ void JacobianBasis::getSignedJacAndGradientsGeneral(int nJacNodes, const fullMat
const fullMatrix<double> &normals,
fullMatrix<double> &JDJ) const
{
-
- switch (bezier->getDim()) {
+ switch (_dim) {
case 0 : {
for (int i = 0; i < nJacNodes; i++) {
@@ -654,9 +652,9 @@ void JacobianBasis::getMetricMinAndGradients(const fullMatrix<double> &nodesXYZ,
gradLambdaJ(l, i + 1 * numMapNodes) = aetaxi * dPhidX + aetaeta * dPhidY;
}
}
-
}
+// Research purpose (to be removed ?)
void JacobianBasis::interpolate(const fullVector<double> &jacobian,
const fullMatrix<double> &uvw,
fullMatrix<double> &result,
@@ -671,7 +669,7 @@ void JacobianBasis::interpolate(const fullVector<double> &jacobian,
else
lag2Bez(jacobian, bez);
- bezier->interpolate(bezM, uvw, result);
+ getBezier()->interpolate(bezM, uvw, result);
}
fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
@@ -683,6 +681,7 @@ fullMatrix<double> JacobianBasis::generateJacMonomialsPyramid(int order)
fullMatrix<double> prox, quad = gmshGenerateMonomialsQuadrangle(2);
prox.setAsProxy(monomials, 0, 2);
prox.setAll(quad);
+ // monomials has been initialized to 0, so third column is ok
return monomials;
}
@@ -789,6 +788,13 @@ fullMatrix<double> JacobianBasis::generateJacPointsPyramid(int order)
return points;
}
+int JacobianBasis::jacobianOrder(int tag)
+{
+ const int parentType = ElementType::ParentTypeFromTag(tag);
+ const int order = ElementType::OrderFromTag(tag);
+ return jacobianOrder(parentType, order);
+}
+
int JacobianBasis::jacobianOrder(int parentType, int order)
{
switch (parentType) {
diff --git a/Numeric/JacobianBasis.h b/Numeric/JacobianBasis.h
index 9701cec635..894b088ce3 100644
--- a/Numeric/JacobianBasis.h
+++ b/Numeric/JacobianBasis.h
@@ -33,6 +33,10 @@ class GradientBasis {
class JacobianBasis {
private:
const GradientBasis *_gradBasis;
+ const bezierBasis *_bezier;
+
+ const int _tag, _dim, _jacOrder;
+
fullMatrix<double> gradShapeMatXFast, gradShapeMatYFast, gradShapeMatZFast;
fullVector<double> primGradShapeBarycenterX, primGradShapeBarycenterY, primGradShapeBarycenterZ;
fullMatrix<double> matrixPrimJac2Jac; // Lifts Lagrange basis of primary Jac. to Lagrange basis of Jac.
@@ -42,19 +46,16 @@ class JacobianBasis {
int numJacNodesFast;
public :
- const bezierBasis *bezier;
-
JacobianBasis(int tag, int jacOrder = -1);
// Get methods
+ inline int getJacOrder() const { return _jacOrder; }
inline int getNumJacNodes() const { return numJacNodes; }
inline int getNumJacNodesFast() const { return numJacNodesFast; }
inline int getNumMapNodes() const { return numMapNodes; }
inline int getNumPrimJacNodes() const { return numPrimJacNodes; }
inline int getNumPrimMapNodes() const { return numPrimMapNodes; }
- inline int getNumDivisions() const { return bezier->getNumDivision(); }
- inline int getNumSubNodes() const { return bezier->subDivisor.size1(); }
- inline int getNumLagCoeff() const { return bezier->getNumLagCoeff(); }
+ const bezierBasis* getBezier() const;
// Jacobian evaluation methods
double getPrimNormals1D(const fullMatrix<double> &nodesXYZ, fullMatrix<double> &result) const;
@@ -99,23 +100,21 @@ class JacobianBasis {
}
//
inline void lag2Bez(const fullVector<double> &jac, fullVector<double> &bez) const {
- bezier->matrixLag2Bez.mult(jac,bez);
+ getBezier()->matrixLag2Bez.mult(jac,bez);
}
inline void lag2Bez(const fullMatrix<double> &jac, fullMatrix<double> &bez) const {
- bezier->matrixLag2Bez.mult(jac,bez);
+ getBezier()->matrixLag2Bez.mult(jac,bez);
}
inline void primJac2Jac(const fullVector<double> &primJac, fullVector<double> &jac) const {
matrixPrimJac2Jac.mult(primJac,jac);
}
- inline void subdivideBezierCoeff(const fullVector<double> &bez, fullVector<double> &result) const {
- bezier->subDivisor.mult(bez,result);
- }
//
void interpolate(const fullVector<double> &jacobian,
const fullMatrix<double> &uvw,
fullMatrix<double> &result, bool areBezier = false) const;
//
+ static int jacobianOrder(int tag);
static int jacobianOrder(int parentType, int order);
static fullMatrix<double> generateJacMonomialsPyramid(int order);
static fullMatrix<double> generateJacPointsPyramid(int order);
diff --git a/Numeric/MetricBasis.cpp b/Numeric/MetricBasis.cpp
index c94a3923ba..093a1d3d86 100644
--- a/Numeric/MetricBasis.cpp
+++ b/Numeric/MetricBasis.cpp
@@ -443,7 +443,7 @@ double MetricBasis::getBoundRmin(MElement *el, MetricData *&md, fullMatrix<doubl
void MetricBasis::_fillInequalities(int metricOrder)
{
- int dimSimplex = _bezier->_dimSimplex;
+ int dimSimplex = _bezier->getDimSimplex();
int dim = _bezier->getDim();
fullMatrix<int> exp(_bezier->_exponents.size1(), _bezier->_exponents.size2());
for (int i = 0; i < _bezier->_exponents.size1(); ++i) {
@@ -533,11 +533,11 @@ void MetricBasis::_fillInequalities(int metricOrder)
}
}
- exp.resize(_jacobian->bezier->_exponents.size1(),
- _jacobian->bezier->_exponents.size2());
- for (int i = 0; i < _jacobian->bezier->_exponents.size1(); ++i) {
- for (int j = 0; j < _jacobian->bezier->_exponents.size2(); ++j) {
- exp(i, j) = static_cast<int>(_jacobian->bezier->_exponents(i, j) + .5);
+ exp.resize(_jacobian->getBezier()->_exponents.size1(),
+ _jacobian->getBezier()->_exponents.size2());
+ for (int i = 0; i < _jacobian->getBezier()->_exponents.size1(); ++i) {
+ for (int j = 0; j < _jacobian->getBezier()->_exponents.size2(); ++j) {
+ exp(i, j) = static_cast<int>(_jacobian->getBezier()->_exponents(i, j) + .5);
}
}
int njac = _jacobian->getNumJacNodes();
@@ -1060,7 +1060,7 @@ double MetricBasis::_subdivideForRmin(
const int numCoeff = md->_metcoeffs->size2();
const int numMetPnts = md->_metcoeffs->size1();
const int numJacPnts = md->_jaccoeffs->size();
- const int numSub = _jacobian->getNumDivisions();
+ const int numSub = _bezier->getNumDivision();
subdomains.push(md);
static unsigned int aa = 200;
@@ -1082,7 +1082,7 @@ double MetricBasis::_subdivideForRmin(
subcoeffs = new fullMatrix<double>(numSub*numMetPnts, numCoeff);
subjac = new fullVector<double>(numSub*numJacPnts);
_bezier->subDivisor.mult(*current->_metcoeffs, *subcoeffs);
- _jacobian->subdivideBezierCoeff(*current->_jaccoeffs, *subjac);
+ _jacobian->getBezier()->subDivisor.mult(*current->_jaccoeffs, *subjac);
int depth = current->_depth;
int num = current->_num;
//Msg::Info("d %d RminBez %g / %g", depth, current->_RminBez, RminLag);
diff --git a/Numeric/bezierBasis.cpp b/Numeric/bezierBasis.cpp
index 0fd6fdb956..edf3952873 100644
--- a/Numeric/bezierBasis.cpp
+++ b/Numeric/bezierBasis.cpp
@@ -13,6 +13,7 @@
#include "BasisFactory.h"
#include "Numeric.h"
+
namespace {
// Sub Control Points
std::vector< fullMatrix<double> > generateSubPointsLine(int order)
@@ -229,8 +230,8 @@ std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
std::vector< fullMatrix<double> > subPoints(4);
fullMatrix<double> prox;
- subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(0);
- subPoints[0].scale(.5/(order+2));
+ subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, 2, 0);
+ subPoints[0].scale(.5/2);
subPoints[1].copy(subPoints[0]);
prox.setAsProxy(subPoints[1], 0, 1);
@@ -250,7 +251,7 @@ std::vector< fullMatrix<double> > generateSubPointsPyr(int order)
std::vector< fullMatrix<double> > subPoints(8);
fullMatrix<double> ref, prox;
- subPoints[0] = JacobianBasis::generateJacMonomialsPyramid(order);
+ subPoints[0] = gmshGenerateMonomialsPyramidGeneral(false, order+2, order);
subPoints[0].scale(.5/(order+2));
subPoints[1].copy(subPoints[0]);
@@ -449,7 +450,7 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
result.resize(uvw.size1(), coeffs.size2());
int dimSimplex;
fullMatrix<double> bezuvw = uvw;
- switch (type) {
+ switch (_type) {
case TYPE_HEX:
if (!bezCoord) {
bezuvw.add(1);
@@ -474,7 +475,7 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
default:
case TYPE_PYR:
- Msg::Error("Bezier interpolation not implemented for type of element %d", type);
+ Msg::Error("Bezier interpolation not implemented for type of element %d", _type);
/*bezuvw[0] = .5 * (1 + uvw[0]);
bezuvw[1] = .5 * (1 + uvw[1]);
bezuvw[2] = uvw[2];
@@ -489,7 +490,7 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
for (int i = 0; i < numCoeff; i++) {
double dd = 1;
double pointCompl = 1.;
- int exponentCompl = order;
+ int exponentCompl = _order;
for (int k = 0; k < dimSimplex; k++) {
dd *= nChoosek(exponentCompl, (int) _exponents(i, k))
* pow(bezuvw(m, k), _exponents(i, k));
@@ -498,10 +499,10 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
}
dd *= pow_int(pointCompl, exponentCompl);
- for (int k = dimSimplex; k < dim; k++) {
- dd *= nChoosek(order, (int) _exponents(i, k))
+ for (int k = dimSimplex; k < _exponents.size2(); k++) {
+ dd *= nChoosek(_order, (int) _exponents(i, k))
* pow_int(bezuvw(m, k), _exponents(i, k))
- * pow_int(1. - bezuvw(m, k), order - _exponents(i, k));
+ * pow_int(1. - bezuvw(m, k), _order - _exponents(i, k));
}
for (int n = 0; n < coeffs.size2(); ++n)
result(m, n) += coeffs(i, n) * dd;
@@ -511,108 +512,99 @@ void bezierBasis::interpolate(const fullMatrix<double> &coeffs,
void bezierBasis::_construct(int parentType, int p)
{
- order = p;
- type = parentType;
+ _order = p;
+ _type = parentType;
std::vector< fullMatrix<double> > subPoints;
if (parentType == TYPE_PYR) {
- dim = 3;
- numLagCoeff = 8;
- _exponents = JacobianBasis::generateJacMonomialsPyramid(order);
+ _numLagCoeff = p == 0 ? 4 : 8;
+ _exponents = gmshGenerateMonomialsPyramidGeneral(false, 2, _order);
- subPoints = generateSubPointsPyr(order);
- numDivisions = static_cast<int>(subPoints.size());
+ subPoints = generateSubPointsPyr(_order);
+ _numDivisions = static_cast<int>(subPoints.size());
fullMatrix<double> bezierPoints;
bezierPoints.resize(_exponents.size1(), _exponents.size2());
- const double p = order + 2;
+ const double ord = _order + 2;
for (int i = 0; i < bezierPoints.size1(); ++i) {
- bezierPoints(i, 2) = _exponents(i, 2) / p;
+ bezierPoints(i, 2) = _exponents(i, 2) / ord;
const double scale = 1. - bezierPoints(i, 2);
- bezierPoints(i, 0) = _exponents(i, 0) / p * scale;
- bezierPoints(i, 1) = _exponents(i, 1) / p * scale;
+ bezierPoints(i, 0) = _exponents(i, 0) / ord * scale;
+ bezierPoints(i, 1) = _exponents(i, 1) / ord * scale;
}
- matrixBez2Lag = generateBez2LagMatrixPyramid(_exponents, bezierPoints, order);
+ matrixBez2Lag = generateBez2LagMatrixPyramid(_exponents, bezierPoints, _order);
matrixBez2Lag.invert(matrixLag2Bez);
- subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez, order);
+ subDivisor = generateSubDivisorPyramid(_exponents, subPoints, matrixLag2Bez, _order);
return;
}
switch (parentType) {
case TYPE_PNT :
- dim = 0;
- numLagCoeff = 1;
+ _numLagCoeff = 1;
_dimSimplex = 0;
_exponents = gmshGenerateMonomialsLine(0);
subPoints.push_back(gmshGeneratePointsLine(0));
break;
case TYPE_LIN : {
- dim = 1;
- numLagCoeff = order == 0 ? 1 : 2;
+ _numLagCoeff = _order == 0 ? 1 : 2;
_dimSimplex = 0;
- _exponents = gmshGenerateMonomialsLine(order);
- subPoints = generateSubPointsLine(order);
+ _exponents = gmshGenerateMonomialsLine(_order);
+ subPoints = generateSubPointsLine(_order);
break;
}
case TYPE_TRI : {
- dim = 2;
- numLagCoeff = order == 0 ? 1 : 3;
+ _numLagCoeff = _order == 0 ? 1 : 3;
_dimSimplex = 2;
- _exponents = gmshGenerateMonomialsTriangle(order);
- subPoints = generateSubPointsTriangle(order);
+ _exponents = gmshGenerateMonomialsTriangle(_order);
+ subPoints = generateSubPointsTriangle(_order);
break;
}
case TYPE_QUA : {
- dim = 2;
- numLagCoeff = order == 0 ? 1 : 4;
+ _numLagCoeff = _order == 0 ? 1 : 4;
_dimSimplex = 0;
- _exponents = gmshGenerateMonomialsQuadrangle(order);
- subPoints = generateSubPointsQuad(order);
+ _exponents = gmshGenerateMonomialsQuadrangle(_order);
+ subPoints = generateSubPointsQuad(_order);
break;
}
case TYPE_TET : {
- dim = 3;
- numLagCoeff = order == 0 ? 1 : 4;
+ _numLagCoeff = _order == 0 ? 1 : 4;
_dimSimplex = 3;
- _exponents = gmshGenerateMonomialsTetrahedron(order);
- subPoints = generateSubPointsTetrahedron(order);
+ _exponents = gmshGenerateMonomialsTetrahedron(_order);
+ subPoints = generateSubPointsTetrahedron(_order);
break;
}
case TYPE_PRI : {
- dim = 3;
- numLagCoeff = order == 0 ? 1 : 6;
+ _numLagCoeff = _order == 0 ? 1 : 6;
_dimSimplex = 2;
- _exponents = gmshGenerateMonomialsPrism(order);
- subPoints = generateSubPointsPrism(order);
+ _exponents = gmshGenerateMonomialsPrism(_order);
+ subPoints = generateSubPointsPrism(_order);
break;
}
case TYPE_HEX : {
- dim = 3;
- numLagCoeff = order == 0 ? 1 : 8;
+ _numLagCoeff = _order == 0 ? 1 : 8;
_dimSimplex = 0;
- _exponents = gmshGenerateMonomialsHexahedron(order);
- subPoints = generateSubPointsHex(order);
+ _exponents = gmshGenerateMonomialsHexahedron(_order);
+ subPoints = generateSubPointsHex(_order);
break;
}
default : {
Msg::Error("Unknown function space of parentType %d : "
"reverting to TET_1", parentType);
- dim = 3;
- order = 0;
- numLagCoeff = 4;
+ _order = 0;
+ _numLagCoeff = 4;
_dimSimplex = 3;
- _exponents = gmshGenerateMonomialsTetrahedron(order);
- subPoints = generateSubPointsTetrahedron(order);
+ _exponents = gmshGenerateMonomialsTetrahedron(_order);
+ subPoints = generateSubPointsTetrahedron(_order);
break;
}
}
- numDivisions = static_cast<int>(subPoints.size());
+ _numDivisions = static_cast<int>(subPoints.size());
fullMatrix<double> bezierPoints = _exponents;
- bezierPoints.scale(1./order);
+ bezierPoints.scale(1./_order);
- matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, order, _dimSimplex);
+ matrixBez2Lag = generateBez2LagMatrix(_exponents, bezierPoints, _order, _dimSimplex);
matrixBez2Lag.invert(matrixLag2Bez);
- subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, order, _dimSimplex);
+ subDivisor = generateSubDivisor(_exponents, subPoints, matrixLag2Bez, _order, _dimSimplex);
}
diff --git a/Numeric/bezierBasis.h b/Numeric/bezierBasis.h
index 257e54be4a..1b0af7216d 100644
--- a/Numeric/bezierBasis.h
+++ b/Numeric/bezierBasis.h
@@ -16,11 +16,11 @@ class MElement;
class bezierBasis {
private :
// the 'numLagCoeff' first exponents are related to 'Lagrangian' values
- int numLagCoeff;
- int dim, type, order;
- int numDivisions;
- public:
- int _dimSimplex;
+ int _numLagCoeff;
+ int _type, _order;
+ int _numDivisions, _dimSimplex;
+
+ friend class MetricBasis;
fullMatrix<double> _exponents;
public :
@@ -37,10 +37,12 @@ class bezierBasis {
}
// get methods
- inline int getDim() const {return dim;}
- inline int getOrder() const {return order;}
- inline int getNumLagCoeff() const {return numLagCoeff;}
- inline int getNumDivision() const {return numDivisions;}
+ inline int getDim() const {return _exponents.size2();}
+ inline int getOrder() const {return _order;}
+ inline int getDimSimplex() const {return _dimSimplex;}
+ inline int getNumLagCoeff() const {return _numLagCoeff;}
+ inline int getNumDivision() const {return _numDivisions;}
+ inline int getNumSubNodes() const {return subDivisor.size1();}
// Interpolation of n functions on N points :
// coeffs(numCoeff, n) and uvw(N, dim)
diff --git a/Numeric/pointsGenerators.cpp b/Numeric/pointsGenerators.cpp
index d8e15c91f2..f060da4ead 100644
--- a/Numeric/pointsGenerators.cpp
+++ b/Numeric/pointsGenerators.cpp
@@ -89,7 +89,7 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
if (order == 0) return points;
for (int i = 0; i < points.size1(); ++i) {
- points(i, 2) = points(i, 2) / order;
+ points(i, 2) = 1 - points(i, 2) / order;
const double duv = -1. + points(i, 2);
points(i, 0) = duv + points(i, 0) * 2. / order;
points(i, 1) = duv + points(i, 1) * 2. / order;
@@ -97,6 +97,23 @@ fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip)
return points;
}
+fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, bool forSerendipPoints)
+{
+ fullMatrix<double> points =
+ gmshGenerateMonomialsPyramidGeneral(pyr, nij, nk, forSerendipPoints);
+
+ if (points.size1() == 1) return points;
+
+ for (int i = 0; i < points.size1(); ++i) {
+ int div = pyr ? nk+nij : std::max(nij, nk);
+ points(i, 2) = (nk - points(i, 2)) / div;
+ const double duv = -1. + points(i, 2);
+ points(i, 0) = duv + points(i, 0) * 2. / div;
+ points(i, 1) = duv + points(i, 1) * 2. / div;
+ }
+ return points;
+}
+
// Monomials Generators
fullMatrix<double> gmshGenerateMonomialsLine(int order, bool serendip)
@@ -668,28 +685,29 @@ fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoint
int nbMonomials = forSerendipPoints ? 5 + (order-1)*8 :
(order+1)*((order+1)+1)*(2*(order+1)+1)/6;
if (forSerendipPoints && !order) nbMonomials = 1;
+
fullMatrix<double> monomials(nbMonomials, 3);
monomials(0, 0) = 0;
monomials(0, 1) = 0;
- monomials(0, 2) = 0;
+ monomials(0, 2) = order;
if (order > 0) {
monomials(1, 0) = order;
monomials(1, 1) = 0;
- monomials(1, 2) = 0;
+ monomials(1, 2) = order;
monomials(2, 0) = order;
monomials(2, 1) = order;
- monomials(2, 2) = 0;
+ monomials(2, 2) = order;
monomials(3, 0) = 0;
monomials(3, 1) = order;
- monomials(3, 2) = 0;
+ monomials(3, 2) = order;
monomials(4, 0) = 0;
monomials(4, 1) = 0;
- monomials(4, 2) = order;
+ monomials(4, 2) = 0;
if (order > 1) {
int index = 5;
@@ -751,7 +769,11 @@ fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoint
if (order > 2) {
fullMatrix<double> inner = gmshGenerateMonomialsPyramid(order - 3);
- inner.add(1);
+ fullMatrix<double> prox;
+ prox.setAsProxy(inner, 0, 2);
+ prox.add(1);
+ prox.setAsProxy(inner, 2, 1);
+ prox.add(2);
monomials.copy(inner, 0, nbMonomials - index, 0, 3, index, 0);
}
}
@@ -760,3 +782,142 @@ fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoint
return monomials;
}
+fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(
+ bool pyr, int nij, int nk, bool forSerendipPoints)
+{
+ if (nij < 0 || nk < 0) {
+ Msg::Fatal("Wrong arguments for pyramid's monomials generation !");
+ }
+ if (!pyr && nk > 0 && nij == 0) {
+ Msg::Error("Wrong argument association for pyramid's monomials generation ! Setting nij to 1");
+ nij = 1;
+ }
+
+ // If monomials for pyramidal space, generate them in gmsh convention.
+ if (forSerendipPoints || (pyr && nij == 0))
+ return gmshGenerateMonomialsPyramid(nk, forSerendipPoints);
+
+ // Otherwise, just put corners at first places,
+ // order of others having no importance.
+
+ int nbMonomials;
+ if (pyr) {
+ nbMonomials = 0;
+ for (int k = 0; k <= nk; ++k) {
+ nbMonomials += (k+nij+1) * (k+nij+1);
+ }
+ }
+ else nbMonomials = (nij+1) * (nij+1) * (nk+1);
+
+ fullMatrix<double> monomials(nbMonomials, 3);
+
+ monomials(0, 0) = 0;
+ monomials(0, 1) = 0;
+ monomials(0, 2) = nk;
+
+ if (nk == 0 && nij == 0) return monomials;
+
+ // Here, nij > 0
+ int nijBase = pyr ? nk+nij : nij;
+
+ // Corners
+ monomials(1, 0) = nijBase;
+ monomials(1, 1) = 0;
+ monomials(1, 2) = nk;
+
+ monomials(2, 0) = nijBase;
+ monomials(2, 1) = nijBase;
+ monomials(2, 2) = nk;
+
+ monomials(3, 0) = 0;
+ monomials(3, 1) = nijBase;
+ monomials(3, 2) = nk;
+
+ if (nk > 0) {
+ monomials(4, 0) = 0;
+ monomials(4, 1) = 0;
+ monomials(4, 2) = 0;
+
+ monomials(5, 0) = nij;
+ monomials(5, 1) = 0;
+ monomials(5, 2) = 0;
+
+ monomials(6, 0) = nij;
+ monomials(6, 1) = nij;
+ monomials(6, 2) = 0;
+
+ monomials(7, 0) = 0;
+ monomials(7, 1) = nij;
+ monomials(7, 2) = 0;
+ }
+
+ int index = 8;
+
+ // Base
+ if (nijBase > 1) {
+ for (int iedge = 0; iedge < 4; ++iedge) {
+ int i0 = iedge;
+ int i1 = (iedge+1) % 4;
+
+ int u0 = (monomials(i1,0)-monomials(i0,0)) / nijBase;
+ int u1 = (monomials(i1,1)-monomials(i0,1)) / nijBase;
+
+ for (int i = 1; i < nijBase; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + i * u0;
+ monomials(index, 1) = monomials(i0, 1) + i * u1;
+ monomials(index, 2) = nk;
+ }
+ }
+
+ for (int i = 1; i < nijBase; ++i) {
+ for (int j = 1; j < nijBase; ++j, ++index) {
+ monomials(index, 0) = i;
+ monomials(index, 1) = j;
+ monomials(index, 2) = nk;
+ }
+ }
+ }
+
+ // Above base
+ if (nk > 0) {
+ // Top
+ if (nij > 1) {
+ for (int iedge = 0; iedge < 4; ++iedge) {
+ int i0 = 4 + iedge;
+ int i1 = 4 + (iedge+1)%4;
+
+ int u0 = (monomials(i1,0)-monomials(i0,0)) / nijBase;
+ int u1 = (monomials(i1,1)-monomials(i0,1)) / nijBase;
+
+ for (int i = 1; i < nijBase; ++i, ++index) {
+ monomials(index, 0) = monomials(i0, 0) + i * u0;
+ monomials(index, 1) = monomials(i0, 1) + i * u1;
+ monomials(index, 2) = 0;
+ }
+ }
+
+ for (int i = 1; i < nijBase; ++i) {
+ for (int j = 1; j < nijBase; ++j, ++index) {
+ monomials(index, 0) = i;
+ monomials(index, 1) = j;
+ monomials(index, 2) = 0;
+ }
+ }
+ }
+
+ // Between bottom & top
+ for (int k = nk-1; k > 0; ++k) {
+ int currentnij = pyr ? k+nij : nij;
+ for (int i = 0; i <= currentnij; ++i) {
+ for (int j = 0; j <= currentnij; ++j, ++index) {
+ monomials(index, 0) = i;
+ monomials(index, 1) = j;
+ monomials(index, 2) = k;
+ }
+ }
+ }
+ }
+
+ return monomials;
+}
+
diff --git a/Numeric/pointsGenerators.h b/Numeric/pointsGenerators.h
index 6546e7dfd6..c3f88f39ab 100644
--- a/Numeric/pointsGenerators.h
+++ b/Numeric/pointsGenerators.h
@@ -29,6 +29,11 @@ fullMatrix<double> gmshGeneratePointsPrism(int order, bool serendip = false);
fullMatrix<double> gmshGeneratePointsPyramid(int order, bool serendip = false);
+// If 'serendip' == true, generate points for serendipity pyramid of order 'nk',
+// else if 'pyr' == true, generate points for space {X^i Y^j Z^k | i,j <= k+'nij', k <= 'nk'},
+// else if 'pyr' == false, generate points for space {X^i Y^j Z^k | i,j <= 'nij', k <= 'nk'}
+fullMatrix<double> gmshGeneratePointsPyramidGeneral(bool pyr, int nij, int nk, bool serendip = false);
+
// Monomial exponents
@@ -45,7 +50,9 @@ fullMatrix<double> gmshGenerateMonomialsPrism(int order, bool forSerendipPoints
fullMatrix<double> gmshGenerateMonomialsPrismSerendipity(int order);
fullMatrix<double> gmshGenerateMonomialsPyramid(int order, bool forSerendipPoints = false);
-//fullMatrix<double> gmshGenerateMonomialsPyramidSerendipity(int order); //TODO
+
+// See explanations for 'gmshGeneratePointsPyramidGeneral'
+fullMatrix<double> gmshGenerateMonomialsPyramidGeneral(bool pyr, int nij, int nk, bool forSerendipPoints = false);
diff --git a/Plugin/AnalyseCurvedMesh.cpp b/Plugin/AnalyseCurvedMesh.cpp
index d3ddac5164..9bf358780e 100644
--- a/Plugin/AnalyseCurvedMesh.cpp
+++ b/Plugin/AnalyseCurvedMesh.cpp
@@ -38,12 +38,12 @@ private:
fullVector<double> _jacBez;
double _minL, _maxL, _minB, _maxB; //Extremum of Jac at corners and of bezier values
int _depthSub;
- const JacobianBasis *_jfs;
+ const bezierBasis *_bfs;
public:
- BezierJacobian(fullVector<double> &, const JacobianBasis *, int depth);
+ BezierJacobian(fullVector<double>&, const bezierBasis*, int depth);
void subDivisions(fullVector<double> &vect) const
- {_jfs->subdivideBezierCoeff(_jacBez, vect);}
+ {_bfs->subDivisor.mult(_jacBez, vect);}
bool boundsOk(double tol, double minL, double maxL) {
return (minL <= 0 || _minB > 0) &&
@@ -294,6 +294,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
+ const bezierBasis *bfs = jfs->getBezier();
_data.reserve(_data.size() + numEl);
@@ -301,7 +302,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el
const int numMapNodes = jfs->getNumMapNodes();
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
- fullVector<double> subJacBez(jfs->getNumSubNodes());
+ fullVector<double> subJacBez(bfs->getNumSubNodes());
for (int k = 0; k < numEl; ++k) {
fullMatrix<double> nodesXYZ(numMapNodes,3);
@@ -309,7 +310,7 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el
jfs->getScaledJacobian(nodesXYZ,jacobian);
jfs->lag2Bez(jacobian, jacBez);
- BezierJacobian *bj = new BezierJacobian(jacBez, jfs, 0);
+ BezierJacobian *bj = new BezierJacobian(jacBez, bfs, 0);
std::vector<BezierJacobian*> heap;
heap.push_back(bj);
@@ -330,9 +331,9 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el
el[k]->getTypeForMSH());
}
- for (int i = 0; i < jfs->getNumDivisions(); i++) {
+ for (int i = 0; i < bfs->getNumDivision(); i++) {
jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
- BezierJacobian *newbj = new BezierJacobian(jacBez, jfs, currentDepth);
+ BezierJacobian *newbj = new BezierJacobian(jacBez, bfs, currentDepth);
minL = std::min(minL, newbj->minL());
maxL = std::max(maxL, newbj->maxL());
@@ -358,9 +359,9 @@ void GMSH_AnalyseCurvedMeshPlugin::_computeMinMaxJandValidity(MElement *const*el
el[k]->getTypeForMSH());
}
- for (int i = 0; i < jfs->getNumDivisions(); i++) {
+ for (int i = 0; i < bfs->getNumDivision(); i++) {
jacBez.setAsProxy(subJacBez, i * numSamplingPt, numSamplingPt);
- BezierJacobian *newbj = new BezierJacobian(jacBez, jfs, currentDepth);
+ BezierJacobian *newbj = new BezierJacobian(jacBez, bfs, currentDepth);
minL = std::min(minL, newbj->minL());
maxL = std::max(maxL, newbj->maxL());
@@ -495,15 +496,15 @@ void GMSH_AnalyseCurvedMeshPlugin::_printStatJacobian()
}
BezierJacobian::BezierJacobian(fullVector<double> &v,
- const JacobianBasis *jfs, int depth)
+ const bezierBasis *bfs, int depth)
{
_jacBez = v;
_depthSub = depth;
- _jfs = jfs;
+ _bfs = bfs;
_minL = _maxL = v(0);
int i = 1;
- for (; i < jfs->getNumLagCoeff(); i++) {
+ for (; i < bfs->getNumLagCoeff(); i++) {
if (_minL > v(i)) _minL = v(i);
if (_maxL < v(i)) _maxL = v(i);
}
diff --git a/Plugin/AnalyseCurvedMesh.h b/Plugin/AnalyseCurvedMesh.h
index 5e9f175cfc..d3b18e445b 100644
--- a/Plugin/AnalyseCurvedMesh.h
+++ b/Plugin/AnalyseCurvedMesh.h
@@ -10,7 +10,6 @@
#include "JacobianBasis.h"
#include "MetricBasis.h"
#include "MElement.h"
-#include "OS.h"
#include <vector>
@@ -81,8 +80,6 @@ public :
void computeMinJ(MElement *const *el, int numEl, double *minJ, bool *straight) {
std::vector<CurvedMeshPluginData> save(_data);
_data.clear();
- double time = Cpu();
- //Rmv add OS.H
_computeMinMaxJandValidity(el, numEl);
if (minJ) {
for (unsigned int i = 0; i < _data.size(); ++i) {
--
GitLab