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

Last change on this file since 8301 was 8301, checked in by Mathieu Morlighem, 14 years ago

moved all grids to nodes
File size: 3.1 KB
Line 
1/*!\file IsInPoly.c
2 * \brief: from a contour, determine which nodes are in domain.
3 */
4
5#include <math.h>
6#include "../../toolkits/toolkits.h"
7#include "./exp.h"
8#include "../shared.h"
9
10#ifdef HAVE_CONFIG_H
11 #include "config.h"
12#else
13#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
14#endif
15
16/*IsInPoly {{{1*/
17int IsInPoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
18
19 int i;
20 double x0,y0;
21 double value;
22 double xmin=xc[0];
23 double xmax=xc[0];
24 double ymin=yc[0];
25 double ymax=yc[0];
26
27 /*Get extrema*/
28 for (i=1;i<numvertices;i++){
29 if(xc[i]<xmin) xmin=xc[i];
30 if(xc[i]>xmax) xmax=xc[i];
31 if(yc[i]<ymin) ymin=yc[i];
32 if(yc[i]>ymax) ymax=yc[i];
33 }
34
35 /*Go through all vertices of the mesh:*/
36 for (i=i0;i<i1;i++){
37
38 //Get current value of value[i] -> do not change it if != 0
39 VecGetValues(in,1,&i,&value);
40 if (value){
41 /*this vertex already is inside one of the contours, continue*/
42 continue;
43 }
44
45 /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
46 x0=x[i]; y0=y[i];
47 if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
48 value=0;
49 }
50 else{
51 value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
52 }
53 VecSetValues(in,1,&i,&value,INSERT_VALUES);
54 }
55 return 1;
56}/*}}}*/
57/*IsOutsidePoly {{{1*/
58int IsOutsidePoly(Vec in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
59
60 int i,j;
61 double x0,y0;
62 double value;
63 double xmin=xc[0];
64 double xmax=xc[0];
65 double ymin=yc[0];
66 double ymax=yc[0];
67
68 /*Get extrema*/
69 for (i=1;i<numvertices;i++){
70 if(xc[i]<xmin) xmin=xc[i];
71 if(xc[i]>xmax) xmax=xc[i];
72 if(yc[i]<ymin) ymin=yc[i];
73 if(yc[i]>ymax) ymax=yc[i];
74 }
75
76 /*Go through all vertices of the mesh:*/
77 for (i=i0;i<i1;i++){
78
79 //Get current value of value[i] -> do not change it if != 0
80 VecGetValues(in,1,&i,&value);
81 if (value){
82 /*this vertex already is inside one of the contours, continue*/
83 continue;
84 }
85
86 /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
87 x0=x[i]; y0=y[i];
88 if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
89 value=1;
90 }
91 else{
92 value=1-pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
93 }
94 VecSetValues(in,1,&i,&value,INSERT_VALUES);
95 }
96 return 1;
97}/*}}}*/
98/*pnpoly{{{1*/
99int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue) {
100 int i, j, c = 0;
101 double n1, n2, normp, scalar;
102
103 /*first test, are they colinear? if yes, is the point in the middle of the segment*/
104 if (edgevalue != 2 ){
105 for (i = 0, j = npol-1; i < npol; j = i++) {
106 n1=pow(yp[i]-yp[j],2.0)+pow(xp[i]-xp[j],2.0);
107 n2=pow(y-yp[j],2.0)+pow(x-xp[j],2.0);
108 normp=pow(n1*n2,0.5);
109 scalar=(yp[i]-yp[j])*(y-yp[j])+(xp[i]-xp[j])*(x-xp[j]);
110
111 if (scalar == normp){
112 if (n2<=n1){
113 c = edgevalue;
114 return c;
115 }
116 }
117 }
118 }
119 /*second test : point is neither on a vertex, nor on a side, where is it ?*/
120 for (i = 0, j = npol-1; i < npol; j = i++) {
121 if ((((yp[i]<=y) && (y<yp[j])) ||
122 ((yp[j]<=y) && (y<yp[i]))) &&
123 (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])){
124 c = !c;
125 }
126 }
127 return c;
128}/*}}}*/
Note: See TracBrowser for help on using the repository browser.