From 53ac59f8ea3a0deb98b8ac7365c55e54bcb61396 Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine <cgeuzaine@ulg.ac.be> Date: Mon, 1 Mar 2010 11:11:32 +0000 Subject: [PATCH] fix msvc compile --- Geo/MQuadrangle.cpp | 135 ++++++++++++++++++++++--------------------- Geo/MTetrahedron.cpp | 20 +++---- Geo/MTriangle.cpp | 7 +-- 3 files changed, 80 insertions(+), 82 deletions(-) diff --git a/Geo/MQuadrangle.cpp b/Geo/MQuadrangle.cpp index eff01bafde..873e3d123e 100644 --- a/Geo/MQuadrangle.cpp +++ b/Geo/MQuadrangle.cpp @@ -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 } diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp index 76b84249da..9c422c1222 100644 --- a/Geo/MTetrahedron.cpp +++ b/Geo/MTetrahedron.cpp @@ -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) diff --git a/Geo/MTriangle.cpp b/Geo/MTriangle.cpp index c2c168cce1..2a060c20e0 100644 --- a/Geo/MTriangle.cpp +++ b/Geo/MTriangle.cpp @@ -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 -- GitLab