Ice Sheet System Model  4.18
Code documentation
Functions
exp.h File Reference
#include <cstring>
#include "../Numerics/recast.h"
#include "../Exceptions/exceptions.h"
#include "../MemOps/MemOps.h"

Go to the source code of this file.

Functions

int IsInPolySerial (double *in, double *xc, double *yc, int numvertices, double *x, double *y, int nods, int edgevalue)
 
int ExpWrite (int nprof, int *profnvertices, double **pprofx, double **pprofy, char *domainname)
 
int pnpoly (int npol, double *xp, double *yp, double x, double y, int edgevalue)
 
template<class doubletype >
int IsInPoly (doubletype *in, double *xc, double *yc, int numvertices, double *x, double *y, int i0, int i1, int edgevalue)
 
template<class doubletype >
int ExpRead (int *pnprof, int **pprofnvertices, doubletype ***ppprofx, doubletype ***ppprofy, bool **pclosed, char *domainname)
 

Function Documentation

◆ IsInPolySerial()

int IsInPolySerial ( double *  in,
double *  xc,
double *  yc,
int  numvertices,
double *  x,
double *  y,
int  nods,
int  edgevalue 
)

Definition at line 45 of file exp.cpp.

45  { /*{{{*/
46 
47  double x0,y0;
48 
49  /*Go through all vertices of the mesh:*/
50  for(int i=0;i<nods;i++){
51  if (in[i]){
52  /*this vertex already is inside one of the contours, continue*/
53  continue;
54  }
55  /*pick up vertex: */
56  x0=x[i];
57  y0=y[i];
58  if (pnpoly(numvertices,xc,yc,x0,y0,edgevalue)){
59  in[i]=1.;
60  }
61  else{
62  in[i]=0.;
63  }
64  }
65 
66  return 1;
67 }

◆ ExpWrite()

int ExpWrite ( int  nprof,
int *  profnvertices,
double **  pprofx,
double **  pprofy,
char *  domainname 
)

Definition at line 13 of file exp.cpp.

13  {/*{{{*/
14 
15  /*I/O: */
16  FILE* fid=NULL;
17 
18  /*open domain outline file for writing: */
19  if((fid=fopen(domainname,"w"))==NULL) _error_("could not open domain file " << domainname);
20 
21  /*Start writing profiles: */
22  for(int counter=0;counter<nprof;counter++){
23 
24  /*Write header: */
25  fprintf(fid,"## Name:%s\n",domainname);
26  fprintf(fid,"## Icon:0\n");
27  fprintf(fid,"# Points Count Value\n");
28  fprintf(fid,"%u %s\n",profnvertices[counter] ,"1.");
29  fprintf(fid,"# X pos Y pos\n");
30 
31  /*Write vertices: */
32  for(int i=0;i<profnvertices[counter];i++){
33  fprintf(fid,"%lf\t%lf\n",pprofx[counter][i],pprofy[counter][i]);
34  }
35 
36  /*Write blank line: */
37  if(counter<nprof-1) fprintf(fid,"\n");
38  }
39 
40  /*close Exp file: */
41  fclose(fid);
42 
43  return 1;
44 }/*}}}*/

◆ pnpoly()

int pnpoly ( int  npol,
double *  xp,
double *  yp,
double  x,
double  y,
int  edgevalue 
)

Definition at line 70 of file exp.cpp.

70  {
71  int i, j, c = 0;
72  double n1, n2, normp, scalar;
73 
74  /*first test, are they colinear? if yes, is the point in the middle of the segment*/
75  if (edgevalue != 2 ){
76  for (i = 0, j = npol-1; i < npol; j = i++) {
77  n1=pow(yp[i]-yp[j],2.0)+pow(xp[i]-xp[j],2.0);
78  n2=pow(y-yp[j],2.0)+pow(x-xp[j],2.0);
79  normp=pow(n1*n2,0.5);
80  scalar=(yp[i]-yp[j])*(y-yp[j])+(xp[i]-xp[j])*(x-xp[j]);
81 
82  if (scalar == normp){
83  if (n2<=n1){
84  c = edgevalue;
85  return c;
86  }
87  }
88  }
89  }
90  /*second test : point is neither on a vertex, nor on a side, where is it ?*/
91  for (i = 0, j = npol-1; i < npol; j = i++) {
92  if ((((yp[i]<=y) && (y<yp[j])) ||
93  ((yp[j]<=y) && (y<yp[i]))) &&
94  (x < (xp[j] - xp[i]) * (y - yp[i]) / (yp[j] - yp[i]) + xp[i])){
95  c = !c;
96  }
97  }
98  return c;
99 }/*}}}*/

◆ IsInPoly()

template<class doubletype >
int IsInPoly ( doubletype *  in,
double *  xc,
double *  yc,
int  numvertices,
double *  x,
double *  y,
int  i0,
int  i1,
int  edgevalue 
)

Definition at line 17 of file exp.h.

17  { /*{{{*/
18 
19  int i;
20  double x0,y0;
21  doubletype 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  value=in[i];
40  if (reCast<bool,doubletype>(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  in[i]=value;
54  }
55  return 1;
56 }/*}}}*/

◆ ExpRead()

template<class doubletype >
int ExpRead ( int *  pnprof,
int **  pprofnvertices,
doubletype ***  ppprofx,
doubletype ***  ppprofy,
bool **  pclosed,
char *  domainname 
)

Definition at line 57 of file exp.h.

57  { /*{{{*/
58 
59  /*indexing: */
60  int i,counter;
61 
62  /*I/O: */
63  FILE *fid = NULL;
64  char chardummy[256];
65  double ddummy;
66 
67  /*output: */
68  int nprof; //number of profiles in the domainname file
69  int *profnvertices = NULL; //array holding the number of vertices for the nprof profiles
70  doubletype **pprofx = NULL; //array of profiles x coordinates
71  doubletype **pprofy = NULL; //array of profiles y coordinates
72  bool *closed = NULL; //array holding closed flags for the nprof profiles
73 
74  /*For each profile: */
75  int n;
76  doubletype *x = NULL;
77  doubletype *y = NULL;
78  bool cl;
79 
80  /*open domain outline file for reading: */
81  if ((fid=fopen(domainname,"r"))==NULL){
82  _error_("could not find file \"" << domainname<<"\". Make sure that the file and path provided exist.");
83  }
84 
85  /*Do a first pass through the domainname file, to figure out how many profiles we need to read: */
86  nprof=1;
87  for(;;){
88  //## Name:filename
89  if(fscanf(fid,"%255s %255s\n",chardummy,chardummy)!=2) _error_("Could not read " << domainname);
90  //## Icon:0
91  if(fscanf(fid,"%255s %255s\n",chardummy,chardummy)!=2) _error_("Could not read " << domainname<<"(Expecting ## Icon:0 and read "<<chardummy<<")");
92  //# Points Count Value
93  if(fscanf(fid,"%255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy)!=4) _error_("Could not read " << domainname);
94  if(fscanf(fid,"%20i %255s\n",&n,chardummy)!=2) _error_("Could not read number of points in "<<domainname);
95  //# X pos Y pos
96  if(fscanf(fid,"%255s %255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy,chardummy)!=5) _error_("Could not read " << domainname);
97  for (i=0;i<n;i++){
98  if(fscanf(fid,"%30lf %30lf\n",&ddummy,&ddummy)!=2){
99  _error_("Could not read coordinate of vertex "<< i <<" of "<<domainname);
100  }
101  }
102  /*check whether we are at the end of the file, otherwise, keep reading next profile:*/
103  if(feof(fid)) break;
104  nprof++;
105  }
106 
107  /*Allocate and initialize all the profiles: */
108  profnvertices = xNew<int>(nprof);
109  pprofx = xNew<doubletype*>(nprof);
110  pprofy = xNew<doubletype*>(nprof);
111  for (i=0;i<nprof;i++){
112  pprofx[i] = NULL;
113  pprofy[i] = NULL;
114  }
115  closed=xNew<bool>(nprof);
116 
117  /*Reaset file pointer to beginning of file: */
118  fseek(fid,0,SEEK_SET);
119 
120  /*Start reading profiles: */
121  for(counter=0;counter<nprof;counter++){
122 
123  /*Skip header: */
124  //## Name:filename
125  if(fscanf(fid,"%255s %255s\n",chardummy,chardummy)!=2) _error_("Could not read " << domainname);
126  //## Icon:0
127  if(fscanf(fid,"%255s %255s\n",chardummy,chardummy)!=2) _error_("Could not read " << domainname);
128  //# Points Count Value
129  if(fscanf(fid,"%255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy)!=4) _error_("Could not read " << domainname);
130 
131  /*Get number of profile vertices: */
132  if(fscanf(fid,"%20i %255s\n",&n,chardummy)!=2) _error_("Could not read number of points in "<<domainname);
133 
134  /*Skip next line: */
135  //# X pos Y pos
136  if(fscanf(fid,"%255s %255s %255s %255s %255s\n",chardummy,chardummy,chardummy,chardummy,chardummy)!=5) _error_("Could not read " << domainname);
137 
138  /*Allocate vertices: */
139  x=xNew<doubletype>(n);
140  y=xNew<doubletype>(n);
141 
142  /*Read vertices: */
143  for (i=0;i<n;i++){
144  if(fscanf(fid,"%30lf %30lf\n",&x[i],&y[i])!=2){
145  _error_("Could not read coordinate of vertex "<<i<<" of "<<domainname);
146  }
147  }
148 
149  /*Now check that we are dealing with open contours: */
150  cl=false;
151  if((x[0]==x[n-1]) && (y[0]==y[n-1])){
152  cl=true;
153  }
154 
155  /*Assign pointers: */
156  profnvertices[counter]=n;
157  pprofx[counter]=x;
158  pprofy[counter]=y;
159  closed[counter]=cl;
160  }
161 
162  /*close domain outline file: */
163  fclose(fid);
164 
165  /*Assign output pointers: */
166  *pnprof=nprof;
167  *pprofnvertices=profnvertices;
168  *ppprofx=pprofx;
169  *ppprofy=pprofy;
170  if(pclosed) *pclosed=closed;
171  else xDelete<bool>(closed);
172  return 1;
173 
174 } /*}}}*/
pnpoly
int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue)
Definition: exp.cpp:70
pnpoly
int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue)
Definition: exp.cpp:70
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19