Ice Sheet System Model  4.18
Code documentation
Functions
InterpFromMeshToGridx.cpp File Reference

"c" core code for interpolating values from a structured grid. More...

#include "./InterpFromMeshToGridx.h"
#include "../../shared/shared.h"

Go to the source code of this file.

Functions

void InterpFromMeshToGridx (double **pgriddata, int *index_mesh, double *x_mesh, double *y_mesh, int nods, int nels, double *data_mesh, int data_length, double *x_grid, double *y_grid, int nlines, int ncols, double default_value)
 

Detailed Description

"c" core code for interpolating values from a structured grid.

Definition in file InterpFromMeshToGridx.cpp.

Function Documentation

◆ InterpFromMeshToGridx()

void InterpFromMeshToGridx ( double **  pgriddata,
int *  index_mesh,
double *  x_mesh,
double *  y_mesh,
int  nods,
int  nels,
double *  data_mesh,
int  data_length,
double *  x_grid,
double *  y_grid,
int  nlines,
int  ncols,
double  default_value 
)

Definition at line 8 of file InterpFromMeshToGridx.cpp.

8  {
9 
10  /*Intermediary*/
11  int i1,i2,j1,j2;
12  int interpolation_type;
13  double area;
14  double area_1,area_2,area_3;
15  double x_tria_min,y_tria_min;
16  double x_tria_max,y_tria_max;
17  double x_grid_min,y_grid_min;
18  double x_grid_max,y_grid_max;
19  double data_value;
20 
21  /*some checks*/
22  if(nels<1 || nods<3 || nlines<2 || ncols<2) _error_("nothing to be done according to the mesh given in input");
23 
24  /*figure out what kind of interpolation is needed*/
25  if(data_length==nods) interpolation_type=1;
26  else if(data_length==nels) interpolation_type=2;
27  else _error_("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
28 
29  /*First, prepare output*/
30  double* griddata=xNew<double>(nlines*ncols);
31  for(int i=0;i<nlines;i++) for(int j=0;j<ncols; j++) griddata[i*ncols+j]=default_value;
32 
33  /*Set debug to "true" if there are lots of elements*/
34  bool debug=(bool)((double)ncols*nlines*nels >= 5.e10);
35 
36  /*figure out if x or y are flipped*/
37  bool xflip = false;
38  bool yflip = false;
39  if(x_grid[1]-x_grid[0]<0) xflip=true;
40  if(y_grid[1]-y_grid[0]<0) yflip=true;
41 
42  /*Get min/max coordinates of the grid*/
43  double xposting = x_grid[1]-x_grid[0];
44  double yposting = y_grid[1]-y_grid[0];
45  if(xflip){
46  x_grid_min=x_grid[ncols-1];
47  x_grid_max=x_grid[0];
48  }
49  else{
50  x_grid_min=x_grid[0];
51  x_grid_max=x_grid[ncols-1];
52  }
53  if(yflip){
54  y_grid_min=y_grid[nlines-1];
55  y_grid_max=y_grid[0];
56  }
57  else{
58  y_grid_min=y_grid[0];
59  y_grid_max=y_grid[nlines-1];
60  }
61 
62  /*Loop over the elements*/
63  for(int n=0;n<nels;n++){
64 
65  /*display current iteration*/
66  if(debug && n%10000==0)
67  _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(n)/double(nels)*100<<"% ");
68 
69  /*Get extrema coordinates of current elements*/
70  x_tria_min=x_mesh[index_mesh[3*n+0]-1]; x_tria_max=x_tria_min;
71  y_tria_min=y_mesh[index_mesh[3*n+0]-1]; y_tria_max=y_tria_min;
72  for(int i=1;i<3;i++){
73  if(x_mesh[index_mesh[3*n+i]-1]<x_tria_min) x_tria_min=x_mesh[index_mesh[3*n+i]-1];
74  if(x_mesh[index_mesh[3*n+i]-1]>x_tria_max) x_tria_max=x_mesh[index_mesh[3*n+i]-1];
75  if(y_mesh[index_mesh[3*n+i]-1]<y_tria_min) y_tria_min=y_mesh[index_mesh[3*n+i]-1];
76  if(y_mesh[index_mesh[3*n+i]-1]>y_tria_max) y_tria_max=y_mesh[index_mesh[3*n+i]-1];
77  }
78 
79  /*if the current triangle is not in the grid, continue*/
80  if( (x_tria_min>x_grid_max) || (x_tria_max<x_grid_min) || (y_tria_min>y_grid_max) || (y_tria_max<y_grid_min) ) continue;
81 
82  /*Get indices i and j that form a square around the current triangle*/
83  if(yflip){
84  i1=max(0, (int)floor((y_tria_max-y_grid_max)/yposting)-1);
85  i2=min(nlines-1,(int)ceil((y_tria_min-y_grid_max)/yposting));
86  }
87  else{
88  i1=max(0, (int)floor((y_tria_min-y_grid_min)/yposting)-1);
89  i2=min(nlines-1,(int)ceil((y_tria_max-y_grid_min)/yposting));
90  }
91  if(xflip){
92  j1=max(0, (int)floor((x_tria_max-x_grid_max)/xposting)-1);
93  j2=min(ncols-1,(int)ceil((x_tria_min-x_grid_max)/xposting));
94  }
95  else{
96  j1=max(0, (int)floor((x_tria_min-x_grid_min)/xposting)-1);
97  j2=min(ncols-1,(int)ceil((x_tria_max-x_grid_min)/xposting));
98  }
99 
100  /*get area of the current element (Jacobian = 2 * area)*/
101  //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
102  area=x_mesh[index_mesh[3*n+1]-1]*y_mesh[index_mesh[3*n+2]-1]-y_mesh[index_mesh[3*n+1]-1]*x_mesh[index_mesh[3*n+2]-1]
103  + x_mesh[index_mesh[3*n+0]-1]*y_mesh[index_mesh[3*n+1]-1]-y_mesh[index_mesh[3*n+0]-1]*x_mesh[index_mesh[3*n+1]-1]
104  + x_mesh[index_mesh[3*n+2]-1]*y_mesh[index_mesh[3*n+0]-1]-y_mesh[index_mesh[3*n+2]-1]*x_mesh[index_mesh[3*n+0]-1];
105 
106  /*Go through x_grid and y_grid and interpolate if necessary*/
107  for(int i=i1;i<=i2;i++){
108 
109  //exit if y not between y_tria_min and y_tria_max
110  if((y_grid[i]>y_tria_max) || (y_grid[i]<y_tria_min)) continue;
111 
112  for(int j=j1;j<=j2; j++){
113 
114  //exit if x not between x_tria_min and x_tria_max
115  if((x_grid[j]>x_tria_max) || (x_grid[j]<x_tria_min)) continue;
116 
117  /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/
118  area_1=((x_grid[j]-x_mesh[index_mesh[3*n+2]-1])*(y_mesh[index_mesh[3*n+1]-1]-y_mesh[index_mesh[3*n+2]-1])
119  - (y_grid[i]-y_mesh[index_mesh[3*n+2]-1])*(x_mesh[index_mesh[3*n+1]-1]-x_mesh[index_mesh[3*n+2]-1]))/area;
120  /*Get second area coordinate =det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/
121  area_2=((x_mesh[index_mesh[3*n+0]-1]-x_mesh[index_mesh[3*n+2]-1])*(y_grid[i]-y_mesh[index_mesh[3*n+2]-1])
122  - (y_mesh[index_mesh[3*n+0]-1]-y_mesh[index_mesh[3*n+2]-1])*(x_grid[j]-x_mesh[index_mesh[3*n+2]-1]))/area;
123  /*Get third area coordinate = 1-area1-area2*/
124  area_3=1.-area_1-area_2;
125 
126  /*is the current point in the current element?*/
127  if(area_1>-1e-12 && area_2>-1e-12 && area_3>-1e-12){
128 
129  /*Yes ! compute the value on the point*/
130  if(interpolation_type==1){
131  /*nodal interpolation*/
132  data_value=area_1*data_mesh[index_mesh[3*n+0]-1]+area_2*data_mesh[index_mesh[3*n+1]-1]+area_3*data_mesh[index_mesh[3*n+2]-1];
133  }
134  else{
135  /*element interpolation*/
136  data_value=data_mesh[n];
137  }
138  if (xIsNan<IssmDouble>(data_value)) data_value=default_value;
139 
140  /*insert value and go to the next point*/
141  griddata[i*ncols+j]=data_value;
142  }
143  }
144  }
145  }
146  if(debug) _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<fixed<<100.<<"% \n");
147 
148  /*Assign output pointers:*/
149  *pgriddata=griddata;
150 }
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
Generated on Thu Jul 2 2020 08:09:23 for Ice Sheet System Model by   doxygen 1.8.19