source: XIOS/trunk/extern/remap/src/intersect.cpp @ 688

#include <cstdlib>
#include <cmath>
#include <cassert>
#include <cstring>
#include <iostream>
#include <fstream>
#include "node.hpp"
#include "elt.hpp"
#include "gridRemap.hpp"
#include "inside.hpp"
#include "polyg.hpp"
#include "intersect.hpp"
#include "intersection_ym.hpp"

namespace sphereRemap {

using namespace std;

/** returns index of edge of a that is shared with b,
    or NOT_FOUND if a and b do not share an edge */
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;
}


/**
If `a` and `b` are neighbours (in the sense that they share an edge)
then this information will be stored in `a` (but not in `b`)
*/
void set_neighbour(Elt& a, const Elt& b)
{
        if (b.id.ind == a.id.ind) return;
        int idx = neighbour_idx(a, b);
        if (idx != NOT_FOUND)
                a.neighbour[idx] = b.id.ind;
}

/** return true if `a` and `b` share an edge */
bool isNeighbour(const Elt& a, const Elt& b)
{
        return neighbour_idx(a, b) != NOT_FOUND;
}

/* computes intersection between elements a and b */

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

        if (!(nseg == nc + nc2 || nseg == 1))
        {


//          cout<<a->x.x<<"  "<<a->x.y<<"  "<<a->x.z<<endl ;
//          cout<<b->x.x<<"  "<<b->x.y<<"  "<<b->x.z<<endl ;
//          cout<<" List of vertex from a and b"<<endl ;
//          for(int i=0;i<na;i++) cout <<"polygonPoints.InsertPoint("<<i<<", "<<a->vertex[i].x<<", "<<a->vertex[i].y<<", "<<a->vertex[i].z<<")"<<endl ;
//          for(int i=0;i<nb;i++) cout <<"polygonPoints.InsertPoint("<<i+6<<", "<<b->vertex[i].x<<", "<<b->vertex[i].y<<", "<<b->vertex[i].z<<")"<<endl ;
//          cout<<"na : "<<na<<"   nb : "<<nb<<endl;
//          cout<<"nc :"<<nc<<"  nc2 :"<<nc2<<"   nseg : "<<nseg<<endl ;
//          abort() ;
//         cout<<"**********************************************"<<endl ;
//         intersect_ym(a,b) ;
        }

//    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]);
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]);
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);
        }

//  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);
/*
          if (  2*fabs(are-ym_area)/(are+ym_area) > 1.1 && ym_area>1e-8)
    {
      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
      intersect_ym(a,b) ;
    }
*/
//              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 )
    {
      cout<<"Big area difference : "<<are<<"  "<<ym_area<<endl ;
      intersect_ym(a,b) ;
    }
*/
//              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)
  {
    if (ym_area>1e-12)
    {
      cout<<"Big area difference : "<<0<<"  "<<ym_area<<endl ;
    }
  }
*/
        delete [] c;
        delete [] c2;
        delete [] xc;
        delete [] xc2;
        delete [] d;
        delete [] d2;
}

}