source: issm/trunk/src/c/shared/Exp/IsInPoly.cpp@ 6965

/*!\file IsInPoly.c
 * \brief:  from a contour, determine which nodes are in  domain.
 */

#include <math.h>
#include "../../toolkits/toolkits.h"
#include "./exp.h"
#include "../shared.h"

#ifdef HAVE_CONFIG_H
        #include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif




int IsInPoly(Vec in,double* xc,double* yc,int numgrids,double* x,double* y,int i0,int i1, int edgevalue){

        int i;
        double x0,y0;
        double value;
        double xmin=xc[0];
        double xmax=xc[0];
        double ymin=yc[0];
        double ymax=yc[0];

        /*Get extrema*/
        for (i=1;i<numgrids;i++){
                if(xc[i]<xmin) xmin=xc[i];
                if(xc[i]>xmax) xmax=xc[i];
                if(yc[i]<ymin) ymin=yc[i];
                if(yc[i]>ymax) ymax=yc[i];
        }

        /*Go through all grids of the mesh:*/
        for (i=i0;i<i1;i++){

                //Get current value of value[i] -> do not change it if != 0
                VecGetValues(in,1,&i,&value);
                if (value){
                        /*this grid already is inside one of the contours, continue*/
                        continue;
                }

                /*pick up grid (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
                x0=x[i]; y0=y[i];
                if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
                        value=0;
                }
                else{
                        value=pnpoly(numgrids,xc,yc,x0,y0,edgevalue);
                }
                VecSetValues(in,1,&i,&value,INSERT_VALUES);
        }
         return 1;
}

int IsOutsidePoly(Vec out,double* xc,double* yc,int numgrids,double* x,double* y,int nods,int edgevalue){

        int i,j;
        double x0,y0;
        double value;

        /*Go through all grids of the mesh:*/
        for (i=MPI_Lowerrow(nods);i<MPI_Upperrow(nods);i++){
                /*pick up grid: */
                x0=x[i];
                y0=y[i];
                if (pnpoly(numgrids,xc,yc,x0,y0,edgevalue)){
                        value=0;
                }
                else{
                        value=1;
                }
                VecSetValues(out,1,&i,&value,INSERT_VALUES);
        }

        return 1;
}

int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue) {
        int i, j, c = 0;
        double n1, n2, normp, scalar;
        
        /*first test, are they colinear? if yes, is the point in the middle of the segment*/
        if (edgevalue != 2 ){
                for (i = 0, j = npol-1; i < npol; j = i++) {
                        n1=pow(yp[i]-yp[j],2.0)+pow(xp[i]-xp[j],2.0);
                        n2=pow(y-yp[j],2.0)+pow(x-xp[j],2.0);
                        normp=pow(n1*n2,0.5);
                        scalar=(yp[i]-yp[j])*(y-yp[j])+(xp[i]-xp[j])*(x-xp[j]);

                        if (scalar == normp){
                                if (n2<=n1){
                                        c = edgevalue;
                                        return c;
                                }
                        }
                }
        }
        /*second test : point is neither on a vertex, nor on a side, where is it ?*/
        for (i = 0, j = npol-1; i < npol; j = i++) {
                if ((((yp[i]<=y) && (y<yp[j])) ||
                                        ((yp[j]<=y) && (y<yp[i]))) &&
                                (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])){
                        c = !c;
                }
        }
        return c;
}