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

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

#include "./InterpFromMeshToMesh3dx.h"
#include "../../shared/shared.h"

Go to the source code of this file.

Functions

int InterpFromMeshToMesh3dx (IssmSeqVec< IssmPDouble > **pdata_prime, double *index_data, double *x_data, double *y_data, double *z_data, int nods_data, int nels_data, double *data, int data_length, double *x_prime, double *y_prime, double *z_prime, int nods_prime, double default_value)
 

Detailed Description

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

Definition in file InterpFromMeshToMesh3dx.cpp.

Function Documentation

◆ InterpFromMeshToMesh3dx()

int InterpFromMeshToMesh3dx ( IssmSeqVec< IssmPDouble > **  pdata_prime,
double *  index_data,
double *  x_data,
double *  y_data,
double *  z_data,
int  nods_data,
int  nels_data,
double *  data,
int  data_length,
double *  x_prime,
double *  y_prime,
double *  z_prime,
int  nods_prime,
double  default_value 
)

Definition at line 8 of file InterpFromMeshToMesh3dx.cpp.

8  {
9 
10  /*Output*/
11  IssmSeqVec<IssmPDouble>* data_prime=NULL;
12 
13  /*Intermediary*/
14  int i,j;
15  int interpolation_type;
16  bool debug;
17  double area;
18  double area_1,area_2,area_3;
19  double zeta,bed,surface;
20  double data_value;
21  double x_prime_min,x_prime_max;
22  double y_prime_min,y_prime_max;
23  double x_tria_min,y_tria_min;
24  double x_tria_max,y_tria_max;
25 
26  /*some checks*/
27  if (nels_data<1 || nods_data<6 || nods_prime==0){
28  _error_("nothing to be done according to the mesh given in input");
29  }
30 
31  /*Set debug to 1 if there are lots of elements*/
32  debug=(bool)((double)nels_data*(double)nods_prime >= pow((double)10,(double)9));
33 
34  /*figure out what kind of interpolation is needed*/
35  if (data_length==nods_data){
36  interpolation_type=1;
37  }
38  else if (data_length==nels_data){
39  interpolation_type=2;
40  }
41  else{
42  _error_("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
43  }
44 
45  /*Get prime mesh extrema coordinates*/
46  x_prime_min=x_prime[0]; x_prime_max=x_prime[0];y_prime_min=y_prime[0]; y_prime_max=y_prime[0];
47  for (i=1;i<nods_prime;i++){
48  if (x_prime[i]<x_prime_min) x_prime_min=x_prime[i];
49  if (x_prime[i]>x_prime_max) x_prime_max=x_prime[i];
50  if (y_prime[i]<y_prime_min) y_prime_min=y_prime[i];
51  if (y_prime[i]>y_prime_max) y_prime_max=y_prime[i];
52  }
53 
54  /*Initialize output*/
55  data_prime=new IssmSeqVec<IssmPDouble>(nods_prime);
56  for (i=0;i<nods_prime;i++) data_prime->SetValue(i,default_value,INS_VAL);
57 
58  /*Loop over the elements*/
59  for (i=0;i<nels_data;i++){
60 
61  /*display current iteration*/
62  if (debug && fmod((double)i,(double)100)==0)
63  _printf_("\r interpolation progress: "<<setw(6)<<setprecision(2)<<double(i)/double(nels_data)*100<<"% ");
64 
65  /*Get extrema coordinates of current elements*/
66  x_tria_min=x_data[(int)index_data[6*i+0]-1]; x_tria_max=x_tria_min;
67  y_tria_min=y_data[(int)index_data[6*i+0]-1]; y_tria_max=y_tria_min;
68  for (j=1;j<3;j++){
69  if(x_data[(int)index_data[6*i+j]-1]<x_tria_min) x_tria_min=x_data[(int)index_data[6*i+j]-1];
70  if(x_data[(int)index_data[6*i+j]-1]>x_tria_max) x_tria_max=x_data[(int)index_data[6*i+j]-1];
71  if(y_data[(int)index_data[6*i+j]-1]<y_tria_min) y_tria_min=y_data[(int)index_data[6*i+j]-1];
72  if(y_data[(int)index_data[6*i+j]-1]>y_tria_max) y_tria_max=y_data[(int)index_data[6*i+j]-1];
73  }
74 
75  /*if there is no point inside the domain, go to next iteration*/
76  if ( x_prime_max < x_tria_min ) continue;
77  if ( x_prime_min > x_tria_max ) continue;
78  if ( y_prime_max < y_tria_min ) continue;
79  if ( y_prime_min > y_tria_max ) continue;
80 
81  /*get area of the current element (Jacobian = 2 * area)*/
82  //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
83  area=x_data[(int)index_data[6*i+1]-1]*y_data[(int)index_data[6*i+2]-1]-y_data[(int)index_data[6*i+1]-1]*x_data[(int)index_data[6*i+2]-1]
84  + x_data[(int)index_data[6*i+0]-1]*y_data[(int)index_data[6*i+1]-1]-y_data[(int)index_data[6*i+0]-1]*x_data[(int)index_data[6*i+1]-1]
85  + x_data[(int)index_data[6*i+2]-1]*y_data[(int)index_data[6*i+0]-1]-y_data[(int)index_data[6*i+2]-1]*x_data[(int)index_data[6*i+0]-1];
86 
87  /*loop over the prime nodes*/
88  for (j=0;j<nods_prime;j++){
89 
90  /*if the current point is not in the triangle, continue*/
91  if ( x_prime[j] < x_tria_min ) continue;
92  if ( x_prime[j] > x_tria_max ) continue;
93  if ( y_prime[j] < y_tria_min ) continue;
94  if ( y_prime[j] > y_tria_max ) continue;
95 
96  /*Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area*/
97  area_1=((x_prime[j]-x_data[(int)index_data[6*i+2]-1])*(y_data[(int)index_data[6*i+1]-1]-y_data[(int)index_data[6*i+2]-1])
98  - (y_prime[j]-y_data[(int)index_data[6*i+2]-1])*(x_data[(int)index_data[6*i+1]-1]-x_data[(int)index_data[6*i+2]-1]))/area;
99  /*Get second area coordinate =det(x1-x3 x-x3 ; y1-y3 y-y3)/area*/
100  area_2=((x_data[(int)index_data[6*i+0]-1]-x_data[(int)index_data[6*i+2]-1])*(y_prime[j]-y_data[(int)index_data[6*i+2]-1])
101  - (y_data[(int)index_data[6*i+0]-1]-y_data[(int)index_data[6*i+2]-1])*(x_prime[j]-x_data[(int)index_data[6*i+2]-1]))/area;
102  /*Get third area coordinate = 1-area1-area2*/
103  area_3=1-area_1-area_2;
104 
105  /*is the current point in the current 2d element?*/
106  if (area_1>=0 && area_2>=0 && area_3>=0){
107 
108  /*compute bottom and top height of the element at this 2d position*/
109  bed =area_1*z_data[(int)index_data[6*i+0]-1]+area_2*z_data[(int)index_data[6*i+1]-1]+area_3*z_data[(int)index_data[6*i+2]-1];
110  surface=area_1*z_data[(int)index_data[6*i+3]-1]+area_2*z_data[(int)index_data[6*i+4]-1]+area_3*z_data[(int)index_data[6*i+5]-1];
111 
112  /*Compute zeta*/
113  zeta=2*(z_prime[j]-bed)/(surface-bed)-1;
114 
115  if (zeta >=-1 && zeta<=1){
116  if (interpolation_type==1){
117  /*nodal interpolation*/
118  data_value=(1-zeta)/2*(area_1*data[(int)index_data[6*i+0]-1]+area_2*data[(int)index_data[6*i+1]-1]+area_3*data[(int)index_data[6*i+2]-1])
119  + (1+zeta)/2*(area_1*data[(int)index_data[6*i+3]-1]+area_2*data[(int)index_data[6*i+4]-1]+area_3*data[(int)index_data[6*i+5]-1]);
120  }
121  else{
122  /*element interpolation*/
123  data_value=data[i];
124  }
125  if (xIsNan<IssmPDouble>(data_value)) data_value=default_value;
126 
127  /*insert value and go to the next point*/
128  data_prime->SetValue(j,data_value,INS_VAL);
129  }
130  }
131  }
132  }
133  if (debug)
134  _printf_("\r interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"% \n");
135 
136  /*Assign output pointers:*/
137  *pdata_prime=data_prime;
138  return 1;
139 }
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
IssmSeqVec< IssmPDouble >
IssmSeqVec::SetValue
void SetValue(int dof, doubletype value, InsMode mode)
Definition: IssmSeqVec.h:115
INS_VAL
@ INS_VAL
Definition: toolkitsenums.h:14
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Generated on Thu Jul 2 2020 08:09:23 for Ice Sheet System Model by   doxygen 1.8.19