Ignore:
Timestamp:
07/20/17 09:18:34 (4 years ago)
Author:
yushan
Message:

test_remap_omp tested on ADA except two fields

Location:
XIOS/dev/branch_openmp/extern/remap/src
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • XIOS/dev/branch_openmp/extern/remap/src/gridRemap.hpp

    r1205 r1220  
    1414Coord readPole(std::istream&); 
    1515 
    16 extern CRemapGrid srcGrid; 
    17 extern CRemapGrid tgtGrid; 
     16 
    1817 
    1918} 
  • XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp

    r1205 r1220  
    1515namespace sphereRemap { 
    1616 
     17extern CRemapGrid srcGrid; 
     18#pragma omp threadprivate(srcGrid) 
     19 
     20extern CRemapGrid tgtGrid; 
     21#pragma omp threadprivate(tgtGrid) 
     22 
    1723using namespace std; 
    1824 
     
    2127int neighbour_idx(const Elt& a, const Elt& b) 
    2228{ 
    23         for (int i = 0; i < a.n; i++) 
    24         { 
    25                 for (int j = 0; j < b.n; j++) 
    26                 { 
    27                         assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS || 
    28                                squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 
    29                         if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 && 
    30                                squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 
    31                         { 
    32                                 return i; 
    33                         } 
    34                 } 
    35         } 
    36         return NOT_FOUND; 
     29  for (int i = 0; i < a.n; i++) 
     30  { 
     31    for (int j = 0; j < b.n; j++) 
     32    { 
     33      assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS || 
     34             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS); 
     35      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 && 
     36             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13) 
     37      { 
     38        return i; 
     39      } 
     40    } 
     41  } 
     42  return NOT_FOUND; 
    3743} 
    3844 
     
    205211bool isNeighbour(Elt& a, const Elt& b) 
    206212{ 
    207         // return neighbour_idx(a, b) != NOT_FOUND; 
     213  // return neighbour_idx(a, b) != NOT_FOUND; 
    208214  return insertNeighbour(a,b,false) ; 
    209215} 
     
    213219void intersect(Elt *a, Elt *b) 
    214220{ 
    215         int na = a->n; /* vertices of a */ 
    216         int nb = b->n; /* vertices of b */ 
    217         Coord *c   = new Coord[na+nb]; 
    218         Coord *c2  = new Coord[na+nb]; 
    219         Coord *xc  = new Coord[na+nb]; 
    220         Coord *xc2 = new Coord[na+nb]; 
    221         Coord gc, gc2; 
    222         double *d = new double[na+nb]; 
    223         double *d2 = new double[na+nb]; 
    224         double are, are2; 
    225         Ipt ipt[NMAX*NMAX]; 
    226         Ipt ipt2[NMAX*NMAX]; 
    227         ptsec(a, b, ipt); 
    228         /* make ipt2 transpose of ipt */ 
    229         for (int ii = 0; ii < na; ii++) 
    230                 for (int jj = 0; jj < nb; jj++) 
    231                         ipt2[jj*na+ii] = ipt[ii*nb+jj]; 
    232         list<Sgm> iscot; 
    233         recense(a, b, ipt, iscot, 0); 
    234         recense(b, a, ipt2, iscot, 1); 
    235  
    236         int nseg = iscot.size(); 
    237         int nc = 0; 
    238         int nc2 = 0; 
    239         while (iscot.size() && nc < 2) 
    240                 nc = assemble(iscot, c, d, xc); 
    241         while (iscot.size() && nc2 < 2) 
    242                 nc2 = assemble(iscot, c2, d2, xc2); 
    243 //      assert(nseg == nc + nc2 || nseg == 1); // unused segment 
     221  int na = a->n; /* vertices of a */ 
     222  int nb = b->n; /* vertices of b */ 
     223  Coord *c   = new Coord[na+nb]; 
     224  Coord *c2  = new Coord[na+nb]; 
     225  Coord *xc  = new Coord[na+nb]; 
     226  Coord *xc2 = new Coord[na+nb]; 
     227  Coord gc, gc2; 
     228  double *d = new double[na+nb]; 
     229  double *d2 = new double[na+nb]; 
     230  double are, are2; 
     231  Ipt ipt[NMAX*NMAX]; 
     232  Ipt ipt2[NMAX*NMAX]; 
     233  ptsec(a, b, ipt); 
     234  /* make ipt2 transpose of ipt */ 
     235  for (int ii = 0; ii < na; ii++) 
     236    for (int jj = 0; jj < nb; jj++) 
     237      ipt2[jj*na+ii] = ipt[ii*nb+jj]; 
     238  list<Sgm> iscot; 
     239  recense(a, b, ipt, iscot, 0); 
     240  recense(b, a, ipt2, iscot, 1); 
     241 
     242  int nseg = iscot.size(); 
     243  int nc = 0; 
     244  int nc2 = 0; 
     245  while (iscot.size() && nc < 2) 
     246    nc = assemble(iscot, c, d, xc); 
     247  while (iscot.size() && nc2 < 2) 
     248    nc2 = assemble(iscot, c2, d2, xc2); 
     249//  assert(nseg == nc + nc2 || nseg == 1); // unused segment 
    244250 
    245251        if (!(nseg == nc + nc2 || nseg == 1)) 
     
    260266 
    261267//    intersect_ym(a,b) ; 
    262         if (nc == 1) nc = 0; 
    263         if (nc2 == 1) nc2 = 0; 
    264         gc = barycentre(xc, nc); 
    265         gc2 = barycentre(xc2, nc2); 
    266         orient(nc, xc, c, d, gc); 
    267  
    268         Coord pole = srcGrid.pole; 
    269         if (pole == ORIGIN) pole = tgtGrid.pole; 
    270         const double MINBASE = 1e-11; 
    271         if (nc == 2) /* nc is the number of vertices of super mesh element */ 
    272         { 
    273                 double base = arcdist(xc[0], xc[1]); 
     268  if (nc == 1) nc = 0; 
     269  if (nc2 == 1) nc2 = 0; 
     270  gc = barycentre(xc, nc); 
     271  gc2 = barycentre(xc2, nc2); 
     272  orient(nc, xc, c, d, gc); 
     273 
     274  Coord pole = srcGrid.pole; 
     275  if (pole == ORIGIN) pole = tgtGrid.pole; 
     276  const double MINBASE = 1e-11; 
     277  if (nc == 2) /* nc is the number of vertices of super mesh element */ 
     278  { 
     279    double base = arcdist(xc[0], xc[1]); 
    274280cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl; 
    275                 gc = midpoint(gc, midpointSC(xc[0], xc[1])); 
    276                 /* intersection area `are` must be zero here unless we have one great and one small circle */ 
    277                 are = alun(base, fabs(scalarprod(xc[0], pole))); 
    278         } 
    279         else 
    280         { 
    281                 are = airbar(nc, xc, c, d, pole, gc); 
    282         } 
    283         if (nc2 == 2) 
    284         { 
    285                 double base = arcdist(xc2[0], xc2[1]); 
     281    gc = midpoint(gc, midpointSC(xc[0], xc[1])); 
     282    /* intersection area `are` must be zero here unless we have one great and one small circle */ 
     283    are = alun(base, fabs(scalarprod(xc[0], pole))); 
     284  } 
     285  else 
     286  { 
     287    are = airbar(nc, xc, c, d, pole, gc); 
     288  } 
     289  if (nc2 == 2) 
     290  { 
     291    double base = arcdist(xc2[0], xc2[1]); 
    286292cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl; 
    287                 assert(base > MINBASE); 
    288                 gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 
    289                 are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 
    290         } 
    291         else 
    292         { 
    293                 are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 
    294         } 
     293    assert(base > MINBASE); 
     294    gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1])); 
     295    are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0 
     296  } 
     297  else 
     298  { 
     299    are2 = airbar(nc2, xc2, c2, d2, pole, gc2); 
     300  } 
    295301 
    296302//  double ym_area=intersect_ym(a,b) ; 
    297303 
    298304  if (nc > 1) 
    299         { 
    300                 /* create one super mesh polygon that src and dest point to */ 
    301                 Polyg *is = new Polyg; 
    302                 is->x = gc; 
    303                 is->area = are; 
    304                 is->id = b->id; 
    305                 is->src_id = b->src_id; 
    306                 is->n = nc; 
    307                 (a->is).push_back(is); 
    308                 (b->is).push_back(is); 
     305  { 
     306    /* create one super mesh polygon that src and dest point to */ 
     307    Polyg *is = new Polyg; 
     308    is->x = gc; 
     309    is->area = are; 
     310    is->id = b->id; 
     311    is->src_id = b->src_id; 
     312    is->n = nc; 
     313    (a->is).push_back(is); 
     314    (b->is).push_back(is); 
    309315/* 
    310           if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 
     316    if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8) 
    311317    { 
    312318      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ; 
     
    314320    } 
    315321*/ 
    316 //              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
    317         } 
    318         if (nc2 > 1) 
    319         { 
    320                 Polyg *is = new Polyg; 
    321                 is->x = gc2; 
    322                 is->area = are2; 
    323                 is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
    324                 is->src_id = b->src_id; 
    325                 is->n = nc2; 
    326                 (a->is).push_back(is); 
    327                 (b->is).push_back(is); 
     322//    cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
     323  } 
     324  if (nc2 > 1) 
     325  { 
     326    Polyg *is = new Polyg; 
     327    is->x = gc2; 
     328    is->area = are2; 
     329    is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */ 
     330    is->src_id = b->src_id; 
     331    is->n = nc2; 
     332    (a->is).push_back(is); 
     333    (b->is).push_back(is); 
    328334/* 
    329     if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 
     335    if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 ) 
    330336    { 
    331337      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ; 
     
    333339    } 
    334340*/ 
    335 //              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
    336         } 
     341//    cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ; 
     342  } 
    337343/* 
    338344  if (nc<=1 && nc2<=1) 
     
    344350  } 
    345351*/ 
    346         delete [] c; 
    347         delete [] c2; 
    348         delete [] xc; 
    349         delete [] xc2; 
    350         delete [] d; 
    351         delete [] d2; 
    352 } 
    353  
    354 } 
     352  delete [] c; 
     353  delete [] c2; 
     354  delete [] xc; 
     355  delete [] xc2; 
     356  delete [] d; 
     357  delete [] d2; 
     358} 
     359 
     360} 
  • XIOS/dev/branch_openmp/extern/remap/src/intersection_ym.cpp

    r1205 r1220  
    1313 
    1414namespace sphereRemap { 
     15 
     16extern CRemapGrid srcGrid; 
     17#pragma omp threadprivate(srcGrid) 
     18 
     19extern CRemapGrid tgtGrid; 
     20#pragma omp threadprivate(tgtGrid) 
     21 
    1522 
    1623using namespace std; 
  • XIOS/dev/branch_openmp/extern/remap/src/mapper.cpp

    r1205 r1220  
    1414 
    1515namespace sphereRemap { 
     16 
     17extern CRemapGrid srcGrid; 
     18#pragma omp threadprivate(srcGrid) 
     19 
     20extern CRemapGrid tgtGrid; 
     21#pragma omp threadprivate(tgtGrid) 
     22 
    1623 
    1724/* A subdivition of an array into N sub-arays 
  • XIOS/dev/branch_openmp/extern/remap/src/parallel_tree.cpp

    r1205 r1220  
    1414 
    1515namespace sphereRemap { 
     16 
     17extern CRemapGrid srcGrid; 
     18#pragma omp threadprivate(srcGrid) 
     19 
     20extern CRemapGrid tgtGrid; 
     21#pragma omp threadprivate(tgtGrid) 
    1622 
    1723static const int assignLevel = 2; 
  • XIOS/dev/branch_openmp/extern/remap/src/polyg.cpp

    r1205 r1220  
    77 
    88#include "polyg.hpp" 
     9#include <stdio.h> 
    910 
    1011namespace sphereRemap { 
     
    4546Coord barycentre(const Coord *x, int n) 
    4647{ 
     48 
    4749        if (n == 0) return ORIGIN; 
    4850        Coord bc = ORIGIN; 
    4951        for (int i = 0; i < n; i++) 
     52        { 
    5053                bc = bc + x[i]; 
     54        }        
    5155        /* both distances can be equal down to roundoff when norm(bc) < mashineepsilon  
    5256           which can occur when weighted with tiny area */ 
    53  
    54         assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0)))); 
    55         //if (squaredist(bc, proj(bc)) > squaredist(bc, proj(bc * (-1.0)))) return proj(bc * (-1.0)); 
     57   
     58         
     59  //assert(squaredist(bc, proj(bc)) <= squaredist(bc, proj(bc * (-1.0))));       
    5660 
    5761        return proj(bc); 
     
    173177                t[1] = x[i]; 
    174178                t[2] = x[ii]; 
    175                 double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 
     179    double sc=scalarprod(crossprod(t[1] - t[0], t[2] - t[0]), t[0]) ; 
    176180                assert(sc >= -1e-10); // Error: tri a l'env (wrong orientation) 
    177181                double area_gc = triarea(t[0], t[1], t[2]); 
     182                if(area_gc<=0) printf("area_gc = %e\n", area_gc); 
    178183                double area_sc_gc_moon = 0; 
    179184                if (d[i]) /* handle small circle case */ 
     
    183188                        char sgl = (mext > 0) ? -1 : 1; 
    184189                        area_sc_gc_moon = sgl * alun(arcdist(t[1], t[2]), fabs(scalarprod(t[1], pole))); 
     190                        //if(area_sc_gc_moon<=0) printf("area_sc_gc_moon = %e\n", area_sc_gc_moon); 
    185191                        gg_exact = gg_exact + sc_gc_moon_normalintegral(t[1], t[2], pole); 
    186192                } 
    187193                area += area_gc + area_sc_gc_moon; /* for "spherical circle segment" sum triangular part (at) and "small moon" and => account for small circle */ 
    188194                g[i] = barycentre(t, 3) * (area_gc + area_sc_gc_moon); 
     195                //printf("g[%d] = (%e,%e,%e) * (%e+%e) = (%e,%e,%e) norm = %e\n", i, barycentre(t, 3).x, barycentre(t, 3).y, barycentre(t, 3).z, area_gc,  area_sc_gc_moon, g[i].x, g[i].y, g[i].z, norm(g[i])); 
    189196        } 
    190197        gg = barycentre(g, N); 
Note: See TracChangeset for help on using the changeset viewer.