Skip to content
Snippets Groups Projects
Commit 53ac59f8 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

fix msvc compile

parent 409bc12b
No related branches found
No related tags found
No related merge requests found
......@@ -189,77 +189,78 @@ double MQuadrangle::angleShapeMeasure()
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 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];
// 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 "
planarQuad_xyz2xy(xp, yp, zp, xn, yn);
// compute the centroid 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;
// 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;
// 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
double R = 1.e22;
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.;
return 0.;
#endif
}
......@@ -49,25 +49,25 @@ double MTetrahedronN::distoShapeMeasure()
double MTetrahedron::getInnerRadius()
{
#if defined(HAVE_MESH)
double r = 0.;
double dist[3];
double face_area = 0.;
double dist[3], face_area = 0.;
double vol = getVolume();
for(int i = 0; i<4; i++){
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;
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]));
}
return 3 * vol / face_area;
#else
r=0.;
return 0.;
#endif
return r;
}
double MTetrahedron::gammaShapeMeasure()
{
#if defined(HAVE_MESH)
......
......@@ -37,16 +37,13 @@ double MTriangle::distoShapeMeasure()
double MTriangle::getInnerRadius()
{
#if defined(HAVE_MESH)
double r = 0.;
double dist[3];
double k = 0.;
double dist[3], k = 0.;
for (int i = 0; i < 3; 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;
return sqrt(k * (k - dist[0]) * (k - dist[1]) * (k - dist[2])) / k;
#else
return 0.;
#endif
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment