From 57253f81793f104c201c0a38756c84ab0254d55a Mon Sep 17 00:00:00 2001
From: Amaury Johnan <amjohnen@gmail.com>
Date: Tue, 10 May 2016 14:45:21 +0000
Subject: [PATCH] common code put in MElement

---
 Geo/MElement.cpp        | 30 +++++++++++++++++++++++++-----
 Geo/MElement.h          |  5 +++--
 Geo/MHexahedron.cpp     |  6 ------
 Geo/MHexahedron.h       |  1 -
 Geo/MLine.cpp           |  6 ------
 Geo/MLine.h             |  1 -
 Geo/MPrism.cpp          |  6 ------
 Geo/MPrism.h            |  1 -
 Geo/MPyramid.cpp        |  6 ------
 Geo/MPyramid.h          |  1 -
 Geo/MQuadrangle.cpp     |  6 ------
 Geo/MQuadrangle.h       |  1 -
 Geo/MTetrahedron.cpp    |  6 ------
 Geo/MTetrahedron.h      |  1 -
 Geo/MTriangle.cpp       |  6 ------
 Geo/MTriangle.h         |  1 -
 Numeric/FuncSpaceData.h |  5 ++---
 17 files changed, 30 insertions(+), 59 deletions(-)

diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index df692d1e28..9c7e895000 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -235,10 +235,10 @@ double MElement::maxDistToStraight() const
   return maxdx;
 }
 
-double MElement::minAnisotropyMeasure()
+double MElement::minAnisotropyMeasure(bool knownValid)
 {
 #if defined(HAVE_MESH)
-  return jacobianBasedQuality::minAnisotropyMeasure(this);
+  return jacobianBasedQuality::minAnisotropyMeasure(this, knownValid);
 #else
   return 0.;
 #endif
@@ -252,7 +252,21 @@ double MElement::specialQuality()
   if (minJ == 0.) return 0;
   if (minJ < 0 && maxJ >= 0) return minJ/maxJ; // accept -inf as an answer
   if (minJ < 0 && maxJ < 0) return -std::numeric_limits<double>::infinity();
-  return jacobianBasedQuality::minAnisotropyMeasure(this);
+  return jacobianBasedQuality::minAnisotropyMeasure(this, true);
+#else
+  return 0;
+#endif
+}
+
+double MElement::specialQuality2()
+{
+#if defined(HAVE_MESH)
+  double minJ, maxJ;
+  jacobianBasedQuality::minMaxJacobianDeterminant(this, minJ, maxJ);
+  if (minJ == 0.) return 0;
+  if (minJ < 0 && maxJ >= 0) return minJ/maxJ; // accept -inf as an answer
+  if (minJ < 0 && maxJ < 0) return -std::numeric_limits<double>::infinity();
+  return jacobianBasedQuality::minScaledJacobian(this, true);
 #else
   return 0;
 #endif
@@ -529,9 +543,9 @@ int MElement::getValidity()
 #if defined(HAVE_MESH)
   double jmin, jmax;
   jacobianBasedQuality::minMaxJacobianDeterminant(this, jmin, jmax);
-  if (jmin > .0 && jmax > .0) return 1; // valid
+  if (jmin > .0) return 1; // valid
   if (jmax >= .0) return 0; // invalid
-  // Here, jmin < 0 and jmax < 0. The element validity is quite indeterminate.
+  // Here, jmax < 0 (and jmin < 0). The element validity is quite indeterminate.
   // It can be valid but with a wrong numbering of the nodes,
   // or it can be invalid, i.e. with nodes that are incorrectly located.
   return -1;
@@ -554,6 +568,12 @@ const nodalBasis* MElement::getFunctionSpace(int order, bool serendip) const
   return tag ? BasisFactory::getNodalBasis(tag) : NULL;
 }
 
+const JacobianBasis* MElement::getJacobianFuncSpace(int order) const
+{
+  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
+  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
+}
+
 static double _computeDeterminantAndRegularize(const MElement *ele, double jac[3][3])
 {
   double dJ = 0;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index ae1b55eb44..6ce208d477 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -216,8 +216,9 @@ class MElement
     signedInvCondNumRange(sICNMin, sICNMax);
     return sICNMin;
   }
-  double minAnisotropyMeasure();
+  double minAnisotropyMeasure(bool knownValid = false);
   double specialQuality();
+  double specialQuality2();
   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);
@@ -268,7 +269,7 @@ class MElement
   virtual const nodalBasis* getFunctionSpace(int order=-1, bool serendip=false) const;
 
   // get the function space for the jacobian of the element
-  virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const { return 0; }
+  virtual const JacobianBasis* getJacobianFuncSpace(int order=-1) const;
 
   // return parametric coordinates (u,v,w) of a vertex
   virtual void getNode(int num, double &u, double &v, double &w) const;
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 630daee66a..9077b1e26a 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -178,12 +178,6 @@ int MHexahedronN::getNumEdgesRep(bool curved)
   return curved ? 12 * CTX::instance()->mesh.numSubEdges : 12;
 }
 
-const JacobianBasis* MHexahedron::getJacobianFuncSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
-}
-
 void _myGetFaceRep(MHexahedron *hex, int num, double *x, double *y, double *z,
                           SVector3 *n, int numSubEdges)
 {
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index 8abbf14c19..7a5fbca9b2 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -59,7 +59,6 @@ class MHexahedron : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual const MVertex *getVertex(int num) const { return _v[num]; }
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual MVertex *getVertexDIFF(int num)
   {
     static const int map[8] = {2, 3, 7, 6, 0, 1, 5, 4};
diff --git a/Geo/MLine.cpp b/Geo/MLine.cpp
index 73e4c79782..4786dc0da2 100644
--- a/Geo/MLine.cpp
+++ b/Geo/MLine.cpp
@@ -13,12 +13,6 @@
 #include "decasteljau.h"
 #include "bezierBasis.h"
 
-const JacobianBasis* MLine::getJacobianFuncSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
-}
-
 void MLine::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
 {
   *npts = getNGQLPts(pOrder);
diff --git a/Geo/MLine.h b/Geo/MLine.h
index 79d4b5d29d..5df4415e15 100644
--- a/Geo/MLine.h
+++ b/Geo/MLine.h
@@ -74,7 +74,6 @@ class MLine : public MElement {
   {
     MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
   }
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual bool isInside(double u, double v, double w) const
   {
     double tol = _isInsideTolerance;
diff --git a/Geo/MPrism.cpp b/Geo/MPrism.cpp
index 12655cc285..488da65857 100644
--- a/Geo/MPrism.cpp
+++ b/Geo/MPrism.cpp
@@ -36,12 +36,6 @@ void MPrism::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *pts = getGQPriPts(pOrder);
 }
 
-const JacobianBasis* MPrism::getJacobianFuncSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
-}
-
 double MPrism::getInnerRadius()
 {
   double dist[3], k = 0.;
diff --git a/Geo/MPrism.h b/Geo/MPrism.h
index 844d147176..c4d4c7cc05 100644
--- a/Geo/MPrism.h
+++ b/Geo/MPrism.h
@@ -128,7 +128,6 @@ class MPrism : public MElement {
     tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
     tmp = _v[3]; _v[3] = _v[4]; _v[4] = tmp;
   }
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual double gammaShapeMeasure();
   virtual int getVolumeSign();
   virtual void getNode(int num, double &u, double &v, double &w) const
diff --git a/Geo/MPyramid.cpp b/Geo/MPyramid.cpp
index b16f096029..5185c65465 100644
--- a/Geo/MPyramid.cpp
+++ b/Geo/MPyramid.cpp
@@ -33,12 +33,6 @@ void MPyramid::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
   *pts = getGQPyrPts(pOrder);
 }
 
-const JacobianBasis* MPyramid::getJacobianFuncSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
-}
-
 MPyramidN::~MPyramidN() {}
 
 int MPyramidN::getNumEdgesRep(bool curved)
diff --git a/Geo/MPyramid.h b/Geo/MPyramid.h
index 85a676b4b7..4e531d8de1 100644
--- a/Geo/MPyramid.h
+++ b/Geo/MPyramid.h
@@ -69,7 +69,6 @@ class MPyramid : public MElement {
   virtual MVertex *getVertex(int num){ return _v[num]; }
   virtual const MVertex *getVertex(int num) const{ return _v[num]; }
   virtual void setVertex(int num,  MVertex *v){ _v[num] = v; }
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual int getNumEdges()const{ return 8; }
   virtual MEdge getEdge(int num) const
   {
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index 1e68dda69f..4216c4d05e 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -16,12 +16,6 @@
 
 #define SQU(a)      ((a)*(a))
 
-const JacobianBasis* MQuadrangle::getJacobianFuncSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
-}
-
 int MQuadrangleN::getNumEdgesRep(bool curved)
 {
   return curved ? 4 * CTX::instance()->mesh.numSubEdges : 4;
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 252507b485..d9abc86b30 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -117,7 +117,6 @@ class MQuadrangle : public MElement {
   virtual const char *getStringForBDF() const { return "CQUAD4"; }
   virtual const char *getStringForDIFF() const { return "ElmB4n2D"; }
   virtual const char *getStringForINP() const { return "CPS4"/*"C2D4"*/; }
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual void getNode(int num, double &u, double &v, double &w) const
   {
     w = 0.;
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 14fd75c1a6..26c57d70a1 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -100,12 +100,6 @@ void MTetrahedron::xyz2uvw(double xyz[3], double uvw[3]) const
   sys3x3(mat, b, uvw, &det);
 }
 
-const JacobianBasis* MTetrahedron::getJacobianFuncSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
-}
-
 int MTetrahedron10::getNumEdgesRep(bool curved){
   return curved ? 6 * CTX::instance()->mesh.numSubEdges : 6;
 }
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index d9ce0ada05..0456dc3959 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -129,7 +129,6 @@ class MTetrahedron : public MElement {
   virtual double getCircumRadius();
   virtual double etaShapeMeasure();
   virtual void xyz2uvw(double xyz[3], double uvw[3]) const;
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual void getNode(int num, double &u, double &v, double &w) const
   {
     switch(num) {
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 3b5638591b..53b3f01729 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -114,12 +114,6 @@ void MTriangle::xyz2uvw(double xyz[3], double uvw[3]) const
   uvw[2] = 0.;
 }
 
-const JacobianBasis* MTriangle::getJacobianFuncSpace(int order) const
-{
-  if (order == -1) return BasisFactory::getJacobianBasis(getTypeForMSH());
-  return BasisFactory::getJacobianBasis(FuncSpaceData(this, order));
-}
-
 int MTriangleN::getNumEdgesRep(bool curved) {
   return curved ? 3 * CTX::instance()->mesh.numSubEdges : 3;
 }
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 3bafc78c3a..08c72ac0ac 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -124,7 +124,6 @@ class MTriangle : public MElement {
   {
     MVertex *tmp = _v[1]; _v[1] = _v[2]; _v[2] = tmp;
   }
-  virtual const JacobianBasis* getJacobianFuncSpace(int o=-1) const;
   virtual void getNode(int num, double &u, double &v, double &w) const
   {
     w = 0.;
diff --git a/Numeric/FuncSpaceData.h b/Numeric/FuncSpaceData.h
index 5f0b889362..b1bb7e57bd 100644
--- a/Numeric/FuncSpaceData.h
+++ b/Numeric/FuncSpaceData.h
@@ -53,9 +53,8 @@ public:
   // Constructors using MElement*
   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);
+  FuncSpaceData(const MElement *el, bool pyr, int nij, int nk,
+                                    const bool *serendip = NULL);
 
   // Constructor using element tag
   FuncSpaceData(int tag, const bool *serendip = NULL);
-- 
GitLab