diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 9afdffeb9757a6876b22e169aa7d29e2f1007fc9..09f188d20309b4e522b195b2dc17fa3bb6b4402d 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -425,7 +425,7 @@ void GFace::computeMeanPlane(const std::vector<SPoint3> &points)
norme(t2);
end:
- res[3] = (xm * res[0] + ym * res[1] + zm * res[2]);
+ res[3] = -(xm * res[0] + ym * res[1] + zm * res[2]);
for(int i = 0; i < 3; i++)
meanPlane.plan[0][i] = t1[i];
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 79a570b27e3f22538ac35a9954ff83e0c32332fe..3a6e92723d815b2c05cb91ba0edc24468b60dfea 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -170,6 +170,11 @@ class MElement
virtual double etaShapeMeasure(){ return 0.; }
virtual double distoShapeMeasure(){ return 1.0; }
virtual double angleShapeMeasure() { return 1.0; }
+
+ // get the radius of the inscribed circle/sphere if it exists,
+ // otherwise get the minimum radius of all the circles/spheres
+ // tangent to the most boundaries of the element.
+ virtual double getInnerRadius(){ return 0.; }
// compute the barycenter
virtual SPoint3 barycenter();
diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp
index fcca6375a34460dcacc01fcc9cf7f1a9d88120ac..eff01bafde6cc39b5f1379894add81dadf58c1ee 100644
--- a/Geo/MQuadrangle.cpp
+++ b/Geo/MQuadrangle.cpp
@@ -7,6 +7,7 @@
#include "GaussLegendre1D.h"
#include "Context.h"
#include "qualityMeasures.h"
+#include <Numeric.h>
#if defined(HAVE_MESH)
#include "qualityMeasures.h"
@@ -185,3 +186,80 @@ double MQuadrangle::angleShapeMeasure()
return 1.;
#endif
}
+
+double MQuadrangle::getInnerRadius()
+{
+ double R = 1.e22;
+
+#if defined(HAVE_MESH)
+ // get the coordinates (x, y, z) of the 4 points defining the Quad
+ double x[4] = {_v[0]->x(), _v[1]->x(), _v[2]->x(), _v[3]->x()};
+ double y[4] = {_v[0]->y(), _v[1]->y(), _v[2]->y(), _v[3]->y()};
+ double z[4] = {_v[0]->z(), _v[1]->z(), _v[2]->z(), _v[3]->z()};
+
+ // get the coefficient (a,b,c,d) of the mean plane - least square!
+ // the plane has for equation " a*x+b*y+c*z+d=0 "
+
+ // compute the centroïd of the 4 points
+ double xm = (x[0]+x[1]+x[2]+x[3])/4;
+ double ym = (y[0]+y[1]+y[2]+y[3])/4;
+ double zm = (z[0]+z[1]+z[2]+z[3])/4;
+
+ // using svd decomposition
+ fullMatrix<double> U(4,3), V(3,3);
+ fullVector<double> sigma(3);
+ for (int i = 0; i < 4; i++) {
+ U(i,0) = x[i]-xm;
+ U(i,1) = y[i]-ym;
+ U(i,2) = z[i]-zm;
+ }
+
+ U.svd(V, sigma);
+ double svd[3];
+ svd[0] = sigma(0);
+ svd[1] = sigma(1);
+ svd[2] = sigma(2);
+ int min;
+ if(fabs(svd[0]) < fabs(svd[1]) && fabs(svd[0]) < fabs(svd[2]))
+ min = 0;
+ else if(fabs(svd[1]) < fabs(svd[0]) && fabs(svd[1]) < fabs(svd[2]))
+ min = 1;
+ else
+ min = 2;
+ double a = V(0, min);
+ double b = V(1, min);
+ double c = V(2, min);
+
+ double d = -(xm * a + ym * b + zm * c);
+
+ double norm = sqrt(a*a+b*b+c*c);
+
+ // projection of the 4 original points on the mean_plane
+
+ double xp[4], yp[4], zp[4];
+
+ for (int i = 0; i < 4; i++) {
+ xp[i] = ((b*b+c*c)*x[i]-a*b*y[i]-a*c*z[i]-d*a)/norm;
+ yp[i] = (-a*b*x[i]+(a*a+c*c)*y[i]-b*c*z[i]-d*b)/norm;
+ zp[i] = (-a*c*x[i]-b*c*y[i]+(a*a+b*b)*z[i]-d*c)/norm;
+ }
+
+ // go from XYZ-plane to XY-plane
+
+ // 4 points, 4 edges => 4 inner radii of circles tangent to (at least) 3 of the four edges!
+ double xn[4], yn[4], r[4];
+
+ planarQuad_xyz2xy(xp, yp, zp, xn, yn);
+
+ // compute for each of the 4 possibilities the incircle radius, keeping the minimum
+ for (int j = 0; j < 4; j++){
+ r[j] = computeInnerRadiusForQuad(xn, yn, j);
+ if(r[j] < R){
+ R = r[j];
+ }
+ }
+ return R;
+#else
+ return 0.;
+#endif
+}
diff --git a/Geo/MQuadrangle.h b/Geo/MQuadrangle.h
index 55e25045447a4b436fcf9a43c1e1b7e68a7db9b0..7ec4a77ff218410c50b21a8bc81db243fbf93df1 100644
--- a/Geo/MQuadrangle.h
+++ b/Geo/MQuadrangle.h
@@ -131,6 +131,10 @@ class MQuadrangle : public MElement {
virtual void getIntegrationPoints(int pOrder, int *npts, IntPt **pts);
virtual double angleShapeMeasure();
virtual double distoShapeMeasure();
+// Computes the minimum inradius of the all the circles tangents to 3 of the 4
+// edges of the quad. If the 4 points of the quad are not planar, we compute
+// the mean plane due to the least-square criterion.
+ virtual double getInnerRadius();
private:
int edges_quad(const int edge, const int vert) const
{
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 10f33374711ce575c96af9953afd7ef54e1b9da9..76b84249da18444b5c11739633d1b0da05e69eab 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -46,6 +46,28 @@ double MTetrahedronN::distoShapeMeasure()
return _disto;
}
+double MTetrahedron::getInnerRadius()
+{
+#if defined(HAVE_MESH)
+ double r = 0.;
+ double dist[3];
+ double face_area = 0.;
+ double vol = getVolume();
+ for(int i = 0; i<4; i++){
+ MFace f = getFace(i);
+ for (int j = 0; j < 3; j++){
+ MEdge e = f.getEdge(j);
+ dist[j] = e.getVertex(0)->distance(e.getVertex(1));
+ }
+ face_area+=0.25*sqrt((dist[0]+dist[1]+dist[2])*(-dist[0]+dist[1]+dist[2])*(dist[0]-dist[1]+dist[2])*(dist[0]+dist[1]-dist[2]));
+ }
+
+r = 3*vol/face_area;
+#else
+ r=0.;
+#endif
+ return r;
+}
double MTetrahedron::gammaShapeMeasure()
{
#if defined(HAVE_MESH)
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index bec81cb6f1c8352b3b9f48b01c1f487ec191ef10..ee917779544b30cb4567dbfb0a0e15f5ee586490 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -128,6 +128,7 @@ class MTetrahedron : public MElement {
virtual int getVolumeSign(){ return (getVolume() >= 0) ? 1 : -1; }
virtual double gammaShapeMeasure();
virtual double distoShapeMeasure();
+ virtual double getInnerRadius(); // radius of sphere inscribed within the Tetrahedron - r=3*Volume/sum(Area_i)
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp
index 4845b4c44b0f05a37d2fcc7be6855cb00f861b94..5e3534c34bcde0020e83f337e9e6f312a3ed04f5 100644
--- a/Geo/MTriangle.cpp
+++ b/Geo/MTriangle.cpp
@@ -34,6 +34,26 @@ double MTriangle::distoShapeMeasure()
#endif
}
+double MTriangle::getInnerRadius()
+{
+#if defined(HAVE_MESH)
+ double r = 0.;
+ int n = getNumEdges();
+ double dist[n];
+ double k = 0.;
+ for (int i=0; i<n; i++)
+ {
+ MEdge e = getEdge(i);
+ dist[i] = e.getVertex(0)->distance(e.getVertex(1));
+ k+=0.5*dist[i];
+ }
+ r=sqrt(k*(k-dist[0])*(k-dist[1])*(k-dist[2]))/k;
+ return r;
+#else
+ return 0.;
+#endif
+}
+
double MTriangle::angleShapeMeasure()
{
#if defined(HAVE_MESH)
diff --git a/Geo/MTriangle.h b/Geo/MTriangle.h
index 30a625539fcbd63f99feac2dba8f451e7d85cfed..ec043e53819eda1fdc9cb4d5a8acfc6e4935096e 100644
--- a/Geo/MTriangle.h
+++ b/Geo/MTriangle.h
@@ -52,6 +52,7 @@ class MTriangle : public MElement {
virtual int getDim() const { return 2; }
virtual double gammaShapeMeasure();
virtual double distoShapeMeasure();
+ virtual double getInnerRadius(); //radius of circle inscribed within the triangle - r=2*Area/sum(Line_i)
virtual double angleShapeMeasure();
virtual int getNumVertices() const { return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }