Changeset 1613


Ignore:
Timestamp:
11/26/18 10:13:26 (21 months ago)
Author:
ymipsl
Message:

interpolation : increase conservation accurency when interpolate integrated flux on the cell. Same methode si now used to compute cell area and intersected cell area.

YM

File:
1 edited

Legend:

Unmodified
Added
Removed
  • XIOS/trunk/extern/remap/src/meshutil.cpp

    r1581 r1613  
    22#include "elt.hpp" 
    33#include "polyg.hpp" 
     4#include "intersection_ym.hpp" 
     5#include "earcut.hpp" 
     6#include <vector> 
    47 
    58namespace sphereRemap { 
    69 
    710using namespace std; 
     11 
     12double 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 
    865 
    966void cptEltGeom(Elt& elt, const Coord &pole) 
     
    1471  elt.area = airbar(elt.n, elt.vertex, elt.edge, elt.d, pole, gg); 
    1572  elt.x = gg; 
    16 } 
     73// overwrite area computation  
     74 
     75  elt.area =  computePolygoneArea(elt, pole) ; 
     76} 
     77 
    1778 
    1879void cptAllEltsGeom(Elt *elt, int N, const Coord &pole) 
Note: See TracChangeset for help on using the changeset viewer.