From 3e536e45652e05004bde1ffd97a196f726c1a99f Mon Sep 17 00:00:00 2001
From: Bruno Seny <bruno.seny@student.uclouvain.be>
Date: Tue, 23 Feb 2010 15:10:44 +0000
Subject: [PATCH] Compute radius of the inscribed circle within an element.
 Done for Mquadrangle, MTriangle, Mtetrahedron

---
 Geo/GFace.cpp        |  2 +-
 Geo/MElement.h       |  5 +++
 Geo/MQuadrangle.cpp  | 78 ++++++++++++++++++++++++++++++++++++++++++++
 Geo/MQuadrangle.h    |  4 +++
 Geo/MTetrahedron.cpp | 22 +++++++++++++
 Geo/MTetrahedron.h   |  1 +
 Geo/MTriangle.cpp    | 20 ++++++++++++
 Geo/MTriangle.h      |  1 +
 8 files changed, 132 insertions(+), 1 deletion(-)

diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index 9afdffeb97..09f188d203 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 79a570b27e..3a6e92723d 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 fcca6375a3..eff01bafde 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 55e2504544..7ec4a77ff2 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 10f3337471..76b84249da 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 bec81cb6f1..ee91777954 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 4845b4c44b..5e3534c34b 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 30a625539f..ec043e5381 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]; }
-- 
GitLab