Skip to content
Snippets Groups Projects
Commit fb548641 authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Changed computation of signed/scaled Jacobian, fixed and completed... 

Changed computation of signed/scaled Jacobian, fixed and completed multi-element versions. Changed "AnalyseCurvedMesh" plugin to work with scaled Jacobian for consistency.
parent ac32ece0
Branches
Tags
No related merge requests found
Show whitespace changes
Inline Side-by-side
Numeric/JacobianBasis.cpp
+ 98
− 101
View file @ fb548641
......@@ -277,27 +277,33 @@ double JacobianBasis::getPrimJac3D(const fullMatrix<double> &nodesXYZ) const
}
namespace {
// Calculate (signed) Jacobian for one element, given the gradients of the shape functions
// and vectors for regularization of line and surface Jacobians in 3D.
inline void computeSignedJacobian(int dim, int numJacNodes, const fullMatrix<double> &gradShapeMatX,
const fullMatrix<double> &gradShapeMatY, const fullMatrix<double> &gradShapeMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
fullVector<double> &jacobian)
// 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>
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
{
const int dim = bezier->getDim();
switch (dim) {
case 0 : {
for (int i = 0; i < numJacNodes; i++) jacobian(i) = 1.;
for (int i = 0; i < nJacNodes; i++) jacobian(i) = 1.;
break;
}
case 1 : {
fullMatrix<double> dxyzdX(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
for (int i = 0; i < numJacNodes; i++) {
fullMatrix<double> normals(2,3);
const double invScale = getPrimNormals1D(nodesXYZ,normals);
if (scaling) {
const double scale = 1./invScale;
normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
}
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);
......@@ -307,29 +313,40 @@ inline void computeSignedJacobian(int dim, int numJacNodes, const fullMatrix<dou
}
case 2 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
gradShapeMatY.mult(nodesXYZ, dxyzdY);
for (int i = 0; i < numJacNodes; i++) {
fullMatrix<double> normal(1,3);
const double invScale = getPrimNormal2D(nodesXYZ,normal);
if (scaling) {
const double scale = 1./invScale;
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
}
fullMatrix<double> dxyzdX(nJacNodes,3), dxyzdY(nJacNodes,3);
gSMatX.mult(nodesXYZ, dxyzdX);
gSMatY.mult(nodesXYZ, dxyzdY);
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);
const double &dxdZ = normal(0,0), &dydZ = normal(0,1), &dzdZ = normal(0,2);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
}
break;
}
case 3 : {
fullMatrix<double> dxyzdX(numJacNodes,3), dxyzdY(numJacNodes,3), dxyzdZ(numJacNodes,3);
gradShapeMatX.mult(nodesXYZ, dxyzdX);
gradShapeMatY.mult(nodesXYZ, dxyzdY);
gradShapeMatZ.mult(nodesXYZ, dxyzdZ);
for (int i = 0; i < numJacNodes; i++) {
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);
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);
jacobian(i) = calcDet3D(dxdX,dydX,dzdX,dxdY,dydY,dzdY,dxdZ,dydZ,dzdZ);
}
if (scaling) {
const double scale = 1./getPrimJac3D(nodesXYZ);
jacobian.scale(scale);
}
break;
}
......@@ -337,88 +354,29 @@ inline void computeSignedJacobian(int dim, int numJacNodes, const fullMatrix<dou
}
} // namespace
// Calculate (signed) Jacobian for one element, with normal vectors to straight element for regularization.
// Evaluation points depend on the given matrices for shape function gradients.
// Calculate signed Jacobian for one element, with normal vectors to straight element for
// regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
const int dim = bezier->getDim();
switch (dim) {
case 1 : {
fullMatrix<double> normals(2,3);
getPrimNormals1D(nodesXYZ,normals);
computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
break;
}
case 2 : {
fullMatrix<double> normal(1,3);
getPrimNormal2D(nodesXYZ,normal);
computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
break;
}
case 0 :
case 3 : {
fullMatrix<double> dum;
computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
break;
getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
}
}
}
// Calculate scaled (signed) Jacobian for one element, with normal vectors to straight element for regularization
// and scaling. Evaluation points depend on the given matrices for shape function gradients.
// Calculate (signed) scaled Jacobian for one element, with normal vectors to straight element
// for regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const
{
const int dim = bezier->getDim();
switch (dim) {
case 1 : {
fullMatrix<double> normals(2,3);
const double scale = 1./getPrimNormals1D(nodesXYZ,normals);
normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normals,jacobian);
break;
}
case 2 : {
fullMatrix<double> normal(1,3);
const double scale = 1./getPrimNormal2D(nodesXYZ,normal);
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,normal,jacobian);
break;
}
case 0 :
case 3 : {
fullMatrix<double> dum;
const double scale = 1./getPrimJac3D(nodesXYZ);
computeSignedJacobian(dim,nJacNodes,gSMatX,gSMatY,gSMatZ,nodesXYZ,dum,jacobian);
jacobian.scale(scale);
break;
}
}
getJacobianGeneral<true>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesXYZ, jacobian);
}
// Calculate (signed) Jacobian for several elements.
// Evaluation points depend on the given matrices for shape function gradients.
// 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
void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
template<bool scaling>
inline void JacobianBasis::getJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
......@@ -438,14 +396,18 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
fullMatrix<double> dxdX(nJacNodes,numEl), dydX(nJacNodes,numEl), dzdX(nJacNodes,numEl);
gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
for (int iEl = 0; iEl < numEl; iEl++) {
fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
for (int i = 0; i < numPrimJacNodes; i++) {
fullMatrix<double> nodesXYZ(numPrimMapNodes,3);
for (int i = 0; i < numPrimMapNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normals(2,3);
getPrimNormals1D(nodesXYZ,normals);
const double invScale = getPrimNormals1D(nodesXYZ,normals);
if (scaling) {
const double scale = 1./invScale;
normals(0,0) *= scale; normals(0,1) *= scale; normals(0,2) *= scale; // Faster to scale 1 normal than afterwards
}
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);
for (int i = 0; i < nJacNodes; i++)
......@@ -463,14 +425,18 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
gSMatX.mult(nodesX, dxdX); gSMatX.mult(nodesY, dydX); gSMatX.mult(nodesZ, dzdX);
gSMatY.mult(nodesX, dxdY); gSMatY.mult(nodesY, dydY); gSMatY.mult(nodesZ, dzdY);
for (int iEl = 0; iEl < numEl; iEl++) {
fullMatrix<double> nodesXYZ(numPrimJacNodes,3);
for (int i = 0; i < numPrimJacNodes; i++) {
fullMatrix<double> nodesXYZ(numPrimMapNodes,3);
for (int i = 0; i < numPrimMapNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
fullMatrix<double> normal(1,3);
getPrimNormal2D(nodesXYZ,normal);
const double invScale = getPrimNormal2D(nodesXYZ,normal);
if (scaling) {
const double scale = 1./invScale;
normal(0,0) *= scale; normal(0,1) *= scale; normal(0,2) *= scale; // Faster to scale normal than afterwards
}
const double &dxdZ = normal(0,0), &dydZ = normal(0,1), &dzdZ = normal(0,2);
for (int i = 0; i < nJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
......@@ -488,11 +454,22 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
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);
for (int iEl = 0; iEl < numEl; iEl++)
for (int iEl = 0; iEl < numEl; iEl++) {
for (int i = 0; i < nJacNodes; i++)
jacobian(i,iEl) = calcDet3D(dxdX(i,iEl),dydX(i,iEl),dzdX(i,iEl),
dxdY(i,iEl),dydY(i,iEl),dzdY(i,iEl),
dxdZ(i,iEl),dydZ(i,iEl),dzdZ(i,iEl));
if (scaling) {
fullMatrix<double> nodesXYZ(numPrimMapNodes,3);
for (int i = 0; i < numPrimMapNodes; i++) {
nodesXYZ(i,0) = nodesX(i,iEl);
nodesXYZ(i,1) = nodesY(i,iEl);
nodesXYZ(i,2) = nodesZ(i,iEl);
}
const double scale = 1./getPrimJac3D(nodesXYZ);
for (int i = 0; i < nJacNodes; i++) jacobian(i,iEl) *= scale;
}
}
break;
}
......@@ -500,6 +477,26 @@ void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<dou
}
// Calculate signed Jacobian for several elements, with normal vectors to straight elements for
// regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
getJacobianGeneral<false>(nJacNodes, gSMatX, gSMatY, gSMatZ, nodesX, nodesY, nodesZ, jacobian);
}
// Calculate (signed) scaled Jacobian for several elements, with normal vectors to straight elements
// for regularization. Evaluation points depend on the given matrices for shape function gradients.
void JacobianBasis::getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const
{
getJacobianGeneral<true>(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,
......
Numeric/JacobianBasis.h
+ 19
− 5
View file @ fb548641
......@@ -86,14 +86,14 @@ class JacobianBasis {
getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesXYZ,jacobian);
}
inline void getScaledJacobian(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodes, _gradBasis->gradShapeMatX,
_gradBasis->gradShapeMatY, _gradBasis->gradShapeMatZ,nodesX,nodesY,nodesZ,jacobian);
}
inline void getSignedJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
}
inline void getSignedJacobianFast(const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const {
getSignedJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,
gradShapeMatZFast,nodesX,nodesY,nodesZ,jacobian);
}
inline void getScaledJacobianFast(const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const {
getScaledJacobianGeneral(numJacNodesFast,gradShapeMatXFast,gradShapeMatYFast,gradShapeMatZFast,nodesXYZ,jacobian);
}
......@@ -122,6 +122,15 @@ class JacobianBasis {
private :
template<bool scaling>
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>
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,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
void getSignedJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
......@@ -132,6 +141,11 @@ class JacobianBasis {
void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, fullVector<double> &jacobian) const;
void getScaledJacobianGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesX, const fullMatrix<double> &nodesY,
const fullMatrix<double> &nodesZ, fullMatrix<double> &jacobian) const;
void getSignedJacAndGradientsGeneral(int nJacNodes, const fullMatrix<double> &gSMatX,
const fullMatrix<double> &gSMatY, const fullMatrix<double> &gSMatZ,
const fullMatrix<double> &nodesXYZ, const fullMatrix<double> &normals,
......
Plugin/AnalyseCurvedMesh.cpp
+ 53
− 90
View file @ fb548641
......@@ -24,15 +24,36 @@
#endif
namespace {
double sum(fullVector<double> &v)
class BezierJacobian;
struct lessMinB {
bool operator()(BezierJacobian*, BezierJacobian*) const;
};
struct lessMaxB {
bool operator()(BezierJacobian*, BezierJacobian*) const;
};
class BezierJacobian
{
double sum = .0;
for (int i = 0; i < v.size(); ++i) {
sum += v(i);
}
return sum;
}
}
private:
fullVector<double> _jacBez;
double _minJ, _maxJ, _minB, _maxB; //Extremum of Jac at corners and of bezier values
int _depthSub;
const JacobianBasis *_jfs;
public:
BezierJacobian(fullVector<double> &, const JacobianBasis *, int depth);
void subDivisions(fullVector<double> &vect) const
{_jfs->subdivideBezierCoeff(_jacBez, vect);}
inline int depth() const {return _depthSub;}
inline double minJ() const {return _minJ;}
inline double maxJ() const {return _maxJ;}
inline double minB() const {return _minB;}
inline double maxB() const {return _maxB;}
};
} // namespace
//#define UNDEF_JAC_TAG -999
//#define _ANALYSECURVEDMESH_BLAS_
......@@ -141,8 +162,6 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
break;
case 2 :
Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
GFace *f = *it;
......@@ -159,7 +178,6 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
GFace *f = *it;
unsigned int numType[3] = {0, 0, 0};
f->getNumMeshElements(numType);
for (int type = 0; type < 3; type++) {
MElement *const *el = f->getStartElementType(type);
checkValidity(el, numType[type], invalids);
......@@ -169,7 +187,6 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(int toDo)
break;
case 1 :
Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
GEdge *e = *it;
unsigned int numElement = e->getNumMeshElements();
......@@ -224,23 +241,20 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
if (numEl < 1)
return;
const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1);
const JacobianBasis *jfs1 = el[0]->getJacobianFuncSpace(1);
if (!jfs || !jfs1) {
const JacobianBasis *jfs = el[0]->getJacobianFuncSpace();
if (!jfs) {
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
const int numSamplingPt = jfs->getNumJacNodes(), numSamplingPt1 = jfs1->getNumJacNodes();
const int numMapNodes = jfs->getNumMapNodes(), numMapNodes1 = jfs1->getNumMapNodes();
const int numSamplingPt = jfs->getNumJacNodes();
const int numMapNodes = jfs->getNumMapNodes();
#ifdef _ANALYSECURVEDMESH_BLAS_
fullMatrix<double> jacobianB(numSamplingPt, numEl);
fullMatrix<double> jacBezB(numSamplingPt, numEl);
fullMatrix<double> jac1B(numSamplingPt1, numEl);
fullVector<double> jacBez, jacobian, jac1;
fullVector<double> jacBez, jacobian;
fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
fullMatrix<double> nodesX1(numMapNodes, numEl), nodesY1(numMapNodes, numEl), nodesZ1(numMapNodes, numEl);
for (int k = 0; k < numEl; ++k)
for (int i = 0; i < numMapNodes; ++i)
{
......@@ -248,20 +262,13 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
nodesX(i,k) = v->x();
nodesY(i,k) = v->y();
nodesZ(i,k) = v->z();
if (i < numMapNodes1) {
nodesX1(i,k) = nodesX(i,k);
nodesY1(i,k) = nodesY(i,k);
nodesZ1(i,k) = nodesZ(i,k);
}
}
jfs->getSignedJacobian(nodesX, nodesY, nodesZ, jacobianB);
jfs1->getSignedJacobian(nodesX1, nodesY1, nodesZ1, jac1B);
jfs->getScaledJacobian(nodesX, nodesY, nodesZ, jacobianB);
jfs->lag2Bez(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
fullVector<double> jac1(numSamplingPt1);
#endif
for (int k = 0; k < numEl; ++k) {
......@@ -269,25 +276,14 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
#ifdef _ANALYSECURVEDMESH_BLAS_
jacBez.setAsProxy(jacBezB, k);
jacobian.setAsProxy(jacobianB, k);
jac1.setAsProxy(jac1B, k);
#else
fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZ1(numMapNodes1,3);
fullMatrix<double> nodesXYZ(numMapNodes,3);
el[k]->getNodesCoord(nodesXYZ);
nodesXYZ1.copy(nodesXYZ,0,numMapNodes1,0,3,0,0);
jfs->getSignedJacobian(nodesXYZ,jacobian);
jfs1->getSignedJacobian(nodesXYZ1,jac1);
jfs->getScaledJacobian(nodesXYZ,jacobian);
#endif
// AmJ : avgJ is not the average Jac for quad, prism or hex
double avgJ = sum(jac1) / jac1.size();
if (avgJ < 0) {
jacBez.scale(-1);
jacobian.scale(-1);
avgJ *= -1;
}
int i;
for (i = 0; i < numSamplingPt && jacobian(i) > _jacBreak * avgJ; ++i);
for (i = 0; i < numSamplingPt && jacobian(i) > _jacBreak; ++i);
if (i < numSamplingPt) {
invalids.push_back(el[k]);
++_numInvalid;
......@@ -304,7 +300,7 @@ void GMSH_AnalyseCurvedMeshPlugin::checkValidity(MElement *const*el,
jfs->lag2Bez(jacobian, jacBez);
#endif
for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak * avgJ; ++i);
for (i = 0; i < jacBez.size() && jacBez(i) > _bezBreak; ++i);
if (i >= jacBez.size()) {
++_numValid;
continue;
......@@ -386,7 +382,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
_numInvalid = 0;
_numValid = 0;
_numUncertain = 0;
_min_Javg = .0, _max_Javg = .0, _avg_Javg = .0;
_min_pJmin = .0, _avg_pJmin = .0;
_min_ratioJ = .0, _avg_ratioJ = .0;
......@@ -406,7 +401,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
break;
case 2 :
Msg::Warning("2D elements must be in a z=cst plane ! If they aren't, results won't be correct.");
for (GModel::fiter it = _m->firstFace(); it != _m->lastFace(); it++) {
GFace *f = *it;
unsigned int numType[3] = {0, 0, 0};
......@@ -421,7 +415,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
break;
case 1 :
Msg::Warning("1D elements must be on a y=cst & z=cst line ! If they aren't, results won't be correct.");
for (GModel::eiter it = _m->firstEdge(); it != _m->lastEdge(); it++) {
GEdge *e = *it;
unsigned int numElement = e->getNumMeshElements();
......@@ -435,7 +428,6 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(std::map<int, std::vector<doubl
Msg::Error("I can't analyse any element.");
return;
}
Msg::Info("Extrema of J_avg : %g, %g (avg: %g)", _min_Javg, _max_Javg, _avg_Javg/_numAnalysedEl);
Msg::Info("Minimum of min(~distortion) : %g", _min_pJmin);
Msg::Info("Average of min(~distortion) : %g", _avg_pJmin / _numAnalysedEl);
Msg::Info("Minimum of min(J) / max(J) : %g", _min_ratioJ);
......@@ -447,24 +439,21 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
if (numEl < 1)
return;
const JacobianBasis *jfs = el[0]->getJacobianFuncSpace(-1);
const JacobianBasis *jfs1 = el[0]->getJacobianFuncSpace(1);
if (!jfs || !jfs1) {
const JacobianBasis *jfs = el[0]->getJacobianFuncSpace();
if (!jfs) {
Msg::Error("Jacobian function space not implemented for type of element %d", el[0]->getNum());
return;
}
const int numSamplingPt = jfs->getNumJacNodes(), numSamplingPt1 = jfs1->getNumJacNodes();
const int numMapNodes = jfs->getNumMapNodes(), numMapNodes1 = jfs1->getNumMapNodes();
const int numSamplingPt = jfs->getNumJacNodes();
const int numMapNodes = jfs->getNumMapNodes();
#ifdef _ANALYSECURVEDMESH_BLAS_
fullMatrix<double> jacobianB(numSamplingPt, numEl);
fullMatrix<double> jacBezB(numSamplingPt, numEl);
fullMatrix<double> jac1B(numSamplingPt1, numEl);
fullVector<double> jacBez, jacobian, jac1;
fullMatrix<double> nodesX(numMapNodes, numEl), nodesY(numMapNodes, numEl), nodesZ(numMapNodes, numEl);
fullMatrix<double> nodesX1(numMapNodes1, numEl), nodesY1(numMapNodes1, numEl), nodesZ1(numMapNodes1, numEl);
for (int k = 0; k < numEl; ++k)
for (int i = 0; i < numMapNodes; ++i)
{
......@@ -472,25 +461,16 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
nodesX(i,k) = v->x();
nodesY(i,k) = v->y();
nodesZ(i,k) = v->z();
if (i < numMapNodes1) {
nodesX1(i,k) = nodesX(i,k);
nodesY1(i,k) = nodesY(i,k);
nodesZ1(i,k) = nodesZ(i,k);
}
}
jfs->getSignedJacobian(nodesX, nodesY, nodesZ, jacobianB);
jfs1->getSignedJacobian(nodesX, nodesY, nodesZ, jac1B);
jfs->getScaledJacobian(nodesX, nodesY, nodesZ, jacobianB);
jfs->lag2Bez(jacobianB, jacBezB);
#else
fullVector<double> jacobian(numSamplingPt);
fullVector<double> jacBez(numSamplingPt);
fullVector<double> jac1(numSamplingPt1);
#endif
fullVector<double> subJacBez(jfs->getNumSubNodes());
_min_Javg = 1.7e308;
_max_Javg = -1.7e308;
_min_pJmin = 1.7e308;
_min_ratioJ = 1.7e308;
......@@ -502,24 +482,13 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
#ifdef _ANALYSECURVEDMESH_BLAS_
jacBez.setAsProxy(jacBezB, k);
jacobian.setAsProxy(jacobianB, k);
jac1.setAsProxy(jac1B, k);
#else
fullMatrix<double> nodesXYZ(numMapNodes,3), nodesXYZ1(numMapNodes1,3);
fullMatrix<double> nodesXYZ(numMapNodes,3);
el[k]->getNodesCoord(nodesXYZ);
nodesXYZ1.copy(nodesXYZ,0,numMapNodes1,0,3,0,0);
jfs->getSignedJacobian(nodesXYZ,jacobian);
jfs1->getSignedJacobian(nodesXYZ1,jac1);
jfs->getScaledJacobian(nodesXYZ,jacobian);
jfs->lag2Bez(jacobian, jacBez);
#endif
// AmJ : avgJ is not the average Jac for quad, prism or hex
double avgJ = sum(jac1) / jac1.size();
if (avgJ < 0) {
jacBez.scale(-1);
jacobian.scale(-1);
avgJ *= -1;
}
double minJ, maxJ = minJ = jacobian(0);
for (int i = 1; i < numSamplingPt; ++i) {
......@@ -527,18 +496,12 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
if (jacobian(i) > maxJ) maxJ = jacobian(i);
}
double minB, maxB = minB = jacBez(0);//, avgJ = .0;
double minB, maxB = minB = jacBez(0);
for (int i = 1; i < numSamplingPt; ++i) {
if (jacBez(i) < minB) minB = jacBez(i);
if (jacBez(i) > maxB) maxB = jacBez(i);
//avgJ += jacBez(i);
}
//avgJ /= numSamplingPt;
_avg_Javg += avgJ;
_min_Javg = std::min(_min_Javg, avgJ);
_max_Javg = std::max(_max_Javg, avgJ);
if (_maxDepth > 1 &&
(minJ - minB > _tol * (std::abs(minJ) + std::abs(minB)) / 2 ||
......@@ -618,18 +581,18 @@ void GMSH_AnalyseCurvedMeshPlugin::computeMinMax(MElement *const*el, int numEl,
delete bj;
}
}
fwrite << minB/avgJ << " " << minB/maxB << "\r";
fwrite << minB << " " << minB/maxB << "\r";
if (data){
if (1-minB <= _tol * minJ && maxB-1 <= _tol * maxB)
(*data)[el[k]->getNum()].push_back(1.);
else if (1-minB/avgJ <= 1e-8)
else if (1-minB <= 1e-8)
(*data)[el[k]->getNum()].push_back(1.);
else
(*data)[el[k]->getNum()].push_back(minB/avgJ);
(*data)[el[k]->getNum()].push_back(minB);
}
_min_pJmin = std::min(_min_pJmin, minB/avgJ);
_avg_pJmin += minB/avgJ;
_min_pJmin = std::min(_min_pJmin, minB);
_avg_pJmin += minB;
_min_ratioJ = std::min(_min_ratioJ, minB/maxB);
_avg_ratioJ += minB/maxB;
}
......
Plugin/AnalyseCurvedMesh.h
+ 0
− 29
View file @ fb548641
......@@ -25,7 +25,6 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
int _numAnalysedEl;
int _numInvalid, _numValid, _numUncertain;
double _min_Javg, _max_Javg, _avg_Javg;
double _min_pJmin, _avg_pJmin;
double _min_ratioJ, _avg_ratioJ;
......@@ -53,32 +52,4 @@ class GMSH_AnalyseCurvedMeshPlugin : public GMSH_PostPlugin
void hideValid_ShowInvalid(std::vector<MElement*> &invalids);
};
class BezierJacobian;
struct lessMinB {
bool operator()(BezierJacobian*, BezierJacobian*) const;
};
struct lessMaxB {
bool operator()(BezierJacobian*, BezierJacobian*) const;
};
class BezierJacobian
{
private:
fullVector<double> _jacBez;
double _minJ, _maxJ, _minB, _maxB; //Extremum of Jac at corners and of bezier values
int _depthSub;
const JacobianBasis *_jfs;
public:
BezierJacobian(fullVector<double> &, const JacobianBasis *, int depth);
void subDivisions(fullVector<double> &vect) const
{_jfs->subdivideBezierCoeff(_jacBez, vect);}
inline int depth() const {return _depthSub;}
inline double minJ() const {return _minJ;}
inline double maxJ() const {return _maxJ;}
inline double minB() const {return _minB;}
inline double maxB() const {return _maxB;}
};
#endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment