Changeset 1613
 Timestamp:
 11/26/18 10:13:26 (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

XIOS/trunk/extern/remap/src/meshutil.cpp
r1581 r1613 2 2 #include "elt.hpp" 3 3 #include "polyg.hpp" 4 #include "intersection_ym.hpp" 5 #include "earcut.hpp" 6 #include <vector> 4 7 5 8 namespace sphereRemap { 6 9 7 10 using namespace std; 11 12 double computePolygoneArea(Elt& a, const Coord &pole) 13 { 14 using N = uint32_t; 15 using Point = array<double, 2>; 16 vector<Point> vect_points; 17 vector< vector<Point> > polyline; 18 19 vector<Coord> dstPolygon ; 20 createGreatCirclePolygon(a, pole, dstPolygon) ; 21 22 int na=dstPolygon.size() ; 23 Coord *a_gno = new Coord[na]; 24 25 Coord OC=barycentre(a.vertex,a.n) ; 26 Coord Oz=OC ; 27 Coord Ox=crossprod(Coord(0,0,1),Oz) ; 28 // choose Ox not too small to avoid rounding error 29 if (norm(Ox)< 0.1) Ox=crossprod(Coord(0,1,0),Oz) ; 30 Ox=Ox*(1./norm(Ox)) ; 31 Coord Oy=crossprod(Oz,Ox) ; 32 double cos_alpha; 33 34 for(int n=0; n<na;n++) 35 { 36 cos_alpha=scalarprod(OC,dstPolygon[n]) ; 37 a_gno[n].x=scalarprod(dstPolygon[n],Ox)/cos_alpha ; 38 a_gno[n].y=scalarprod(dstPolygon[n],Oy)/cos_alpha ; 39 a_gno[n].z=scalarprod(dstPolygon[n],Oz)/cos_alpha ; // must be equal to 1 40 41 vect_points.push_back( array<double, 2>() ); 42 vect_points[n][0] = a_gno[n].x; 43 vect_points[n][1] = a_gno[n].y; 44 45 } 46 47 polyline.push_back(vect_points); 48 vector<N> indices_a_gno = mapbox::earcut<N>(polyline); 49 50 double area_a_gno=0 ; 51 for(int i=0;i<indices_a_gno.size()/3;++i) 52 { 53 Coord x0 = Ox * polyline[0][indices_a_gno[3*i]][0] + Oy* polyline[0][indices_a_gno[3*i]][1] + Oz ; 54 Coord x1 = Ox * polyline[0][indices_a_gno[3*i+1]][0] + Oy* polyline[0][indices_a_gno[3*i+1]][1] + Oz ; 55 Coord x2 = Ox * polyline[0][indices_a_gno[3*i+2]][0] + Oy* polyline[0][indices_a_gno[3*i+2]][1] + Oz ; 56 area_a_gno+=triarea(x0 * (1./norm(x0)),x1* (1./norm(x1)), x2* (1./norm(x2))) ; 57 } 58 59 vect_points.clear(); 60 polyline.clear(); 61 indices_a_gno.clear(); 62 return area_a_gno ; 63 } 64 8 65 9 66 void cptEltGeom(Elt& elt, const Coord &pole) … … 14 71 elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg); 15 72 elt.x = gg; 16 } 73 // overwrite area computation 74 75 elt.area = computePolygoneArea(elt, pole) ; 76 } 77 17 78 18 79 void cptAllEltsGeom(Elt *elt, int N, const Coord &pole)
Note: See TracChangeset
for help on using the changeset viewer.