Commit fa80be93 by Kirill Terekhov

Another fix to geometric algorithms

parent e639c57d
 ... ... @@ -1099,6 +1099,8 @@ namespace INMOST vec_diff(v1,cxf,l1,mdim); vec_diff(v2,cxf,l2,mdim); vec_cross_product(l1,l2,nt); for(int q = 0; q < 3; ++q) nt[q] *= 0.5; ss = vec_dot_product(n0,nt,3); if( ss ) ss /= fabs(ss); at = sqrt(vec_dot_product(nt,nt,3))*ss; ... ... @@ -1123,10 +1125,11 @@ namespace INMOST rev = CreatePrivateMarker(); FacesOrientation(faces,rev); } real vol = 0, a, at; real vol = 0, a, at, volp; real x[3] = {0,0,0}, n[3] = {0,0,0}, n0[3] = {0,0,0}, s, ss; real l1[3] = {0,0,0}, l2[3] = {0,0,0}; real nt[3] = {0,0,0}, cxf[3] = {0,0,0}; real vc[3] = {0,0,0}; //me.Centroid(cx); //if( print ) std::cout << "cx: " << cx[0] << " " << cx[1] << " " << cx[2] << std::endl; for(unsigned j = 0; j < faces.size(); j++) ... ... @@ -1141,13 +1144,8 @@ namespace INMOST n[0] = n[1] = n[2] = 0; n0[0] = n0[1] = n0[2] = 0; cxf[0] = cxf[1] = cxf[2] = 0; vc[0] = vc[1] = vc[2] = 0; a = 0; //real_array v0 = nodes[0].Coords(); //real_array v1 = nodes[1].Coords(); //real_array v2 = nodes[2].Coords(); //vec_diff(v1,v0,l1,mdim); //vec_diff(v2,v0,l2,mdim); //vec_cross_product(l1,l2,n0); real_array /*v0 = nodes[0].Coords(),*/ v1, v2; faces[j].Centroid(cxf); for(int i = 0; i < (int)nodes.size(); i++) ... ... @@ -1157,11 +1155,9 @@ namespace INMOST vec_diff(v1,cxf,l1,mdim); vec_diff(v2,cxf,l2,mdim); vec_cross_product(l1,l2,nt); //if( print ) std::cout << "nt: " << nt[0] << " " << nt[1] << " " << nt[2] << std::endl; for(int q = 0; q < 3; ++q) n0[q] += nt[q]*0.5; } //if( print ) std::cout << "n0: " << n0[0] << " " << n0[1] << " " << n0[2] << std::endl; for(int k = 0; k < (int)nodes.size(); k++) { v1 = nodes[k].Coords(); ... ... @@ -1174,7 +1170,6 @@ namespace INMOST ss = vec_dot_product(n0,nt,3); if( ss ) ss /= fabs(ss); at = sqrt(vec_dot_product(nt,nt,3))*ss; //if( print ) std::cout << "nt: " << nt[0] << " " << nt[1] << " " << nt[2] << " ss " << ss << " at " << at << std::endl; //same as faces[j].Normal(n) for(int q = 0; q < 3; ++q) n[q] += nt[q]; ... ... @@ -1183,33 +1178,22 @@ namespace INMOST x[q] += at*(/*(v0[q]-cx[q])+*/(v1[q]-cxf[q])+(v2[q]-cxf[q]))/3.0; a += at; } //if( print ) std::cout << "n: " << n[0] << " " << n[1] << " " << n[2] << std::endl; //if( print ) std::cout << "x: " << x[0] << " " << x[1] << " " << x[2] << " a " << a << std::endl; if( a ) { //~ std::cout << "Correction: " << x[0]/a << " " << x[1]/a << " " << x[2]/a << " area " << a << std::endl; for(int q = 0; q < 3; ++q) x[q] = x[q]/a + cxf[q]; //~ std::cout << "Centroid: " << cxf[0] << " " << cxf[1] << " " << cxf[2] << std::endl; //~ std::cout << "Barycenter: " << x[0] << " " << x[1] << " " << x[2] << std::endl; } else for(int q = 0; q < 3; ++q) x[q] = cxf[q]; vol += s*vec_dot_product(x,n,3); volp = vec_dot_product(x,n,3) / 3.0; vol += s*volp; } //if( print ) std::cout << "vol: " << vol << std::endl; if( ornt ) { if( vol < 0.0 ) vol = -vol; faces.RemPrivateMarker(rev); ReleasePrivateMarker(rev); } *ret = vol/3.0; //if( *ret < 0 || *ret != *ret ) //{ // std::cout << "volume is " << *ret << std::endl; // print = true; // goto redo; //} *ret = vol; break; } } ... ... @@ -1270,12 +1254,6 @@ namespace INMOST *ret = 0; real nt[3] = {0,0,0}, l1[3] = {0,0,0}, l2[3] = {0,0,0}; real c[3] = {0,0,0}, n0[3] = {0,0,0}, ss, cxf[3] = {0,0,0}; //real_array v0 = nodes[0].Coords(); //real_array v1 = nodes[1].Coords(); //real_array v2 = nodes[2].Coords(); //vec_diff(v1,v0,l1,mdim); //vec_diff(v2,v0,l2,mdim); //vec_cross_product(l1,l2,n0); real_array v0 = nodes[0].Coords(), v1, v2; Element(this,e).Centroid(cxf); for(int i = 0; i < (int)nodes.size(); i++) ... ... @@ -1296,9 +1274,11 @@ namespace INMOST vec_diff(v1,cxf,l1,mdim); vec_diff(v2,cxf,l2,mdim); vec_cross_product(l1,l2,nt); for(int q = 0; q < 3; ++q) nt[q] *= 0.5; ss = vec_dot_product(n0,nt,3); if( ss ) ss /= fabs(ss); at = sqrt(vec_dot_product(nt,nt,3))*0.5*ss; at = sqrt(vec_dot_product(nt,nt,3))*ss; for(int q = 0; q < mdim; ++q) c[q] += at*(/*(v0[q]-cxf[q])+*/(v1[q]-cxf[q])+(v2[q]-cxf[q]))/3.0; a += at; ... ... @@ -1309,7 +1289,6 @@ namespace INMOST ret[q] = c[q]/a+cxf[q]; } else for(int q = 0; q < mdim; ++q) ret[q] = cxf[q]; } //std::cout << ret[0] << " " << ret[1] << " " << ret[2] << std::endl; } else if( edim == 3 ) { ... ... @@ -1322,12 +1301,13 @@ namespace INMOST rev = CreatePrivateMarker(); FacesOrientation(faces,rev); } real vol = 0, a, at; real vol = 0, a, at, volp; real x[3] = {0,0,0}, nt[3] = {0,0,0}, s; real c[3] = {0,0,0}, n[3] = {0,0,0}; real n0[3] = {0,0,0}, ss; real l1[3] = {0,0,0}, l2[3] = {0,0,0}; real cx[3] = {0,0,0}, cxf[3] = {0,0,0}; real vc[3] = {0,0,0}; me.Centroid(cx); for(unsigned j = 0; j < faces.size(); j++) { ... ... @@ -1337,16 +1317,11 @@ namespace INMOST s = faces[j].GetPrivateMarker(rev) ? -1.0 : 1.0; else s = faces[j].FaceOrientedOutside(me) ? 1.0 : -1.0; //real_array v0 = nodes[0].Coords(); //real_array v1 = nodes[1].Coords(); //real_array v2 = nodes[2].Coords(); //vec_diff(v1,v0,l1,mdim); //vec_diff(v2,v0,l2,mdim); //vec_cross_product(l1,l2,n0); n0[0] = n0[1] = n0[2] = 0; x[0] = x[1] = x[2] = 0; n[0] = n[1] = n[2] = 0; cxf[0] = cxf[1] = cxf[2] = 0; vc[0] = vc[1] = vc[2] = 0; a = 0; real_array /*v0 = nodes[0].Coords(),*/ v1, v2; faces[j].Centroid(cxf); ... ... @@ -1379,16 +1354,17 @@ namespace INMOST for(int q = 0; q < 3; ++q) x[q] += at*(/*(v0[q]-cxf[q])+*/(v1[q]-cxf[q])+(v2[q]-cxf[q]))/3.0; a += at; //second-order midpoint formula for(int q = 0; q < 3; ++q) c[q] += s*nt[q]*(pow((/*v0*/cxf[q]-cx[q])+(v1[q]-cx[q]),2)+pow((/*v0*/cxf[q]-cx[q])+(v2[q]-cx[q]),2)+pow((v1[q]-cx[q])+(v2[q]-cx[q]),2))/24.0; } if( a ) { for(int q = 0; q < 3; ++q) x[q] = x[q]/a + cxf[q]; } else for(int q = 0; q < 3; ++q) x[q] = cxf[q]; vol += s*vec_dot_product(x,n,3); vec_diff(x,cx,vc,3); volp = vec_dot_product(vc,n,3) / 3.0; for(int q = 0; q < 3; ++q) c[q] += 0.75*(x[q]-cx[q])*s*volp; // x - base barycenter, cx - apex vol += s*volp; } if( ornt ) { ... ... @@ -1401,7 +1377,6 @@ namespace INMOST faces.RemPrivateMarker(rev); ReleasePrivateMarker(rev); } vol /= 3.0; if( vol ) { for(int q = 0; q < mdim; ++q) ... ...
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!