Changeset 1328 for XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp

Timestamp:

11/15/17 12:14:34 (7 years ago)

Author:

yushan

Message:

dev_omp

File:

• XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp (modified) (8 diffs)

XIOS/dev/branch_openmp/extern/remap/src/intersect.cpp

@@ -15 +15 @@
 namespace sphereRemap {
 
-extern CRemapGrid srcGrid;
-#pragma omp threadprivate(srcGrid)
-
-extern CRemapGrid tgtGrid;
-#pragma omp threadprivate(tgtGrid)
-
 using namespace std;
 
@@ -27 +21 @@
 int neighbour_idx(const Elt& a, const Elt& b)
 {
-  for (int i = 0; i < a.n; i++)
-  {
-    for (int j = 0; j < b.n; j++)
-    {
-      assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS ||
-             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS);
-      if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 &&
-             squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13)
-      {
-        return i;
-      }
-    }
-  }
-  return NOT_FOUND;
+        for (int i = 0; i < a.n; i++)
+        {
+                for (int j = 0; j < b.n; j++)
+                {
+                        assert(squaredist(a.vertex[ i       ], b.vertex[ j       ]) > EPS*EPS ||
+                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+1)%b.n]) > EPS*EPS);
+                        if (   squaredist(a.vertex[ i       ], b.vertex[ j           ]) < 1e-13*1e-13 &&
+                               squaredist(a.vertex[(i+1)%a.n], b.vertex[(j+b.n-1)%b.n]) < 1e-13*1e-13)
+                        {
+                                return i;
+                        }
+                }
+        }
+        return NOT_FOUND;
 }
 
@@ -211 +205 @@
 bool isNeighbour(Elt& a, const Elt& b)
 {
-  // return neighbour_idx(a, b) != NOT_FOUND;
+        // return neighbour_idx(a, b) != NOT_FOUND;
   return insertNeighbour(a,b,false) ;
 }
@@ -219 +213 @@
 void intersect(Elt *a, Elt *b)
 {
-  int na = a->n; /* vertices of a */
-  int nb = b->n; /* vertices of b */
-  Coord *c   = new Coord[na+nb];
-  Coord *c2  = new Coord[na+nb];
-  Coord *xc  = new Coord[na+nb];
-  Coord *xc2 = new Coord[na+nb];
-  Coord gc, gc2;
-  double *d = new double[na+nb];
-  double *d2 = new double[na+nb];
-  double are, are2;
-  Ipt ipt[NMAX*NMAX];
-  Ipt ipt2[NMAX*NMAX];
-  ptsec(a, b, ipt);
-  /* make ipt2 transpose of ipt */
-  for (int ii = 0; ii < na; ii++)
-    for (int jj = 0; jj < nb; jj++)
-      ipt2[jj*na+ii] = ipt[ii*nb+jj];
-  list<Sgm> iscot;
-  recense(a, b, ipt, iscot, 0);
-  recense(b, a, ipt2, iscot, 1);
-
-  int nseg = iscot.size();
-  int nc = 0;
-  int nc2 = 0;
-  while (iscot.size() && nc < 2)
-    nc = assemble(iscot, c, d, xc);
-  while (iscot.size() && nc2 < 2)
-    nc2 = assemble(iscot, c2, d2, xc2);
-//  assert(nseg == nc + nc2 || nseg == 1); // unused segment
+        int na = a->n; /* vertices of a */
+        int nb = b->n; /* vertices of b */
+        Coord *c   = new Coord[na+nb];
+        Coord *c2  = new Coord[na+nb];
+        Coord *xc  = new Coord[na+nb];
+        Coord *xc2 = new Coord[na+nb];
+        Coord gc, gc2;
+        double *d = new double[na+nb];
+        double *d2 = new double[na+nb];
+        double are, are2;
+        Ipt ipt[NMAX*NMAX];
+        Ipt ipt2[NMAX*NMAX];
+        ptsec(a, b, ipt);
+        /* make ipt2 transpose of ipt */
+        for (int ii = 0; ii < na; ii++)
+                for (int jj = 0; jj < nb; jj++)
+                        ipt2[jj*na+ii] = ipt[ii*nb+jj];
+        list<Sgm> iscot;
+        recense(a, b, ipt, iscot, 0);
+        recense(b, a, ipt2, iscot, 1);
+
+        int nseg = iscot.size();
+        int nc = 0;
+        int nc2 = 0;
+        while (iscot.size() && nc < 2)
+                nc = assemble(iscot, c, d, xc);
+        while (iscot.size() && nc2 < 2)
+                nc2 = assemble(iscot, c2, d2, xc2);
+//      assert(nseg == nc + nc2 || nseg == 1); // unused segment
 
         if (!(nseg == nc + nc2 || nseg == 1))
@@ -266 +260 @@
 
 //    intersect_ym(a,b) ;
-  if (nc == 1) nc = 0;
-  if (nc2 == 1) nc2 = 0;
-  gc = barycentre(xc, nc);
-  gc2 = barycentre(xc2, nc2);
-  orient(nc, xc, c, d, gc);
-
-  Coord pole = srcGrid.pole;
-  if (pole == ORIGIN) pole = tgtGrid.pole;
-  const double MINBASE = 1e-11;
-  if (nc == 2) /* nc is the number of vertices of super mesh element */
-  {
-    double base = arcdist(xc[0], xc[1]);
+        if (nc == 1) nc = 0;
+        if (nc2 == 1) nc2 = 0;
+        gc = barycentre(xc, nc);
+        gc2 = barycentre(xc2, nc2);
+        orient(nc, xc, c, d, gc);
+
+        Coord pole = srcGrid.pole;
+        if (pole == ORIGIN) pole = tgtGrid.pole;
+        const double MINBASE = 1e-11;
+        if (nc == 2) /* nc is the number of vertices of super mesh element */
+        {
+                double base = arcdist(xc[0], xc[1]);
 cerr << "DID ARRIVE " << base << xc[0] << xc[1] << endl;
-    gc = midpoint(gc, midpointSC(xc[0], xc[1]));
-    /* intersection area `are` must be zero here unless we have one great and one small circle */
-    are = alun(base, fabs(scalarprod(xc[0], pole)));
-  }
-  else
-  {
-    are = airbar(nc, xc, c, d, pole, gc);
-  }
-  if (nc2 == 2)
-  {
-    double base = arcdist(xc2[0], xc2[1]);
+                gc = midpoint(gc, midpointSC(xc[0], xc[1]));
+                /* intersection area `are` must be zero here unless we have one great and one small circle */
+                are = alun(base, fabs(scalarprod(xc[0], pole)));
+        }
+        else
+        {
+                are = airbar(nc, xc, c, d, pole, gc);
+        }
+        if (nc2 == 2)
+        {
+                double base = arcdist(xc2[0], xc2[1]);
 cerr << "DID ARRIVE " << base << xc2[0] << xc2[1] << endl;
-    assert(base > MINBASE);
-    gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1]));
-    are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0
-  }
-  else
-  {
-    are2 = airbar(nc2, xc2, c2, d2, pole, gc2);
-  }
+                assert(base > MINBASE);
+                gc2 = midpoint(gc2, midpointSC(xc2[0], xc2[1]));
+                are2 = alun(base, fabs(scalarprod(xc2[0], pole))); // 0
+        }
+        else
+        {
+                are2 = airbar(nc2, xc2, c2, d2, pole, gc2);
+        }
 
 //  double ym_area=intersect_ym(a,b) ;
 
   if (nc > 1)
-  {
-    /* create one super mesh polygon that src and dest point to */
-    Polyg *is = new Polyg;
-    is->x = gc;
-    is->area = are;
-    is->id = b->id;
-    is->src_id = b->src_id;
-    is->n = nc;
-    (a->is).push_back(is);
-    (b->is).push_back(is);
+        {
+                /* create one super mesh polygon that src and dest point to */
+                Polyg *is = new Polyg;
+                is->x = gc;
+                is->area = are;
+                is->id = b->id;
+                is->src_id = b->src_id;
+                is->n = nc;
+                (a->is).push_back(is);
+                (b->is).push_back(is);
 /*
-    if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
+          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
     {
       cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
@@ -320 +314 @@
     }
 */
-//    cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
-  }
-  if (nc2 > 1)
-  {
-    Polyg *is = new Polyg;
-    is->x = gc2;
-    is->area = are2;
-    is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
-    is->src_id = b->src_id;
-    is->n = nc2;
-    (a->is).push_back(is);
-    (b->is).push_back(is);
+//              cout<<"intersection : "<<are<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
+        }
+        if (nc2 > 1)
+        {
+                Polyg *is = new Polyg;
+                is->x = gc2;
+                is->area = are2;
+                is->id = b->id; /* intersection holds id of corresponding source element (see Elt class definition for details about id) */
+                is->src_id = b->src_id;
+                is->n = nc2;
+                (a->is).push_back(is);
+                (b->is).push_back(is);
 /*
-    if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 )
+    if (        2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8 )
     {
       cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
@@ -339 +333 @@
     }
 */
-//    cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
-  }
+//              cout<<"intersection : "<<are2<<" "<< ym_area<<"  diff : "<<fabs(are-ym_area)<<"  ratio : "<<fabs(are-ym_area)/(0.5*(are+ym_area))<<endl ;
+        }
 /*
   if (nc<=1 && nc2<=1)
@@ -350 +344 @@
   }
 */
-  delete [] c;
-  delete [] c2;
-  delete [] xc;
-  delete [] xc2;
-  delete [] d;
-  delete [] d2;
-}
-
-}
+        delete [] c;
+        delete [] c2;
+        delete [] xc;
+        delete [] xc2;
+        delete [] d;
+        delete [] d2;
+}
+
+}