InterpFromMesh2dx.h File Reference

: header file for Data interpolation routines. More...

#include "../../classes/classes.h"
#include "../../toolkits/toolkits.h"

Go to the source code of this file.

Data Structures
struct	InterpFromMesh2dxThreadStruct

Functions
int	InterpFromMesh2dx (IssmSeqVec< IssmPDouble > **pdata_prime, int *index_data, double *x_data, double *y_data, int nods_data, int nels_data, double *data, int data_length, double *x_prime, double *y_prime, int nods_prime, double *default_values, int num_default_values, Contour< IssmPDouble > **contours, int numcontours)
void *	InterpFromMesh2dxt (void *vInterpFromMesh2dxThreadStruct)

Detailed Description

: header file for Data interpolation routines.

Definition in file InterpFromMesh2dx.h.

Function Documentation

InterpFromMesh2dx()

int InterpFromMesh2dx	(	IssmSeqVec< IssmPDouble > **	pdata_prime,
		int *	index_data,
		double *	x_data,
		double *	y_data,
		int	nods_data,
		int	nels_data,
		double *	data,
		int	data_length,
		double *	x_prime,
		double *	y_prime,
		int	nods_prime,
		double *	default_values,
		int	num_default_values,
		Contour< IssmPDouble > **	contours,
		int	numcontours
	)

Definition at line 11 of file InterpFromMesh2dx.cpp.

                                                                                                       {
 
   /*Output*/
   IssmSeqVec<IssmPDouble>* data_prime=NULL;
 
   /*Intermediary*/
   int    i;
   int    interpolation_type;
   bool   debug;
   double xmin,xmax;
   double ymin,ymax;
 
   /*contours: */
   double *incontour     = NULL;
 
   /*some checks*/
   if (nels_data<1 || nods_data<3 || nods_prime==0){
      _error_("nothing to be done according to the mesh given in input");
   }
 
   /*Set debug to 1 if there are lots of elements*/
   debug=(bool)((double)nels_data*(double)nods_prime >= pow((double)10,(double)9));
 
   /*figure out what kind of interpolation is needed*/
   if (data_length==nods_data){
      interpolation_type=1;
   }
   else if (data_length==nels_data){
      interpolation_type=2;
   }
   else{
      _error_("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
   }
 
   if((numcontours) && (interpolation_type==2)){
      _error_("element interpolation_type with contours not supported yet!");
   }
 
   /*Get prime mesh extrema coordinates*/
   xmin=x_prime[0]; xmax=x_prime[0];ymin=y_prime[0]; ymax=y_prime[0];
   for (i=1;i<nods_prime;i++){
      if (x_prime[i]<xmin) xmin=x_prime[i];
      if (x_prime[i]>xmax) xmax=x_prime[i];
      if (y_prime[i]<ymin) ymin=y_prime[i];
      if (y_prime[i]>ymax) ymax=y_prime[i];
   }
 
   /*Initialize output*/
   data_prime=new IssmSeqVec<IssmPDouble>(nods_prime);
   if(num_default_values){
      if(num_default_values==1)for (i=0;i<nods_prime;i++) data_prime->SetValue(i,default_values[0],INS_VAL);
      else for (i=0;i<nods_prime;i++) data_prime->SetValue(i,default_values[i],INS_VAL);
   }
 
   /*Build indices of contour: */
   if(numcontours){
      ContourToNodesx( &incontour,x_prime,y_prime,nods_prime,contours,numcontours,1);
   }
   else{
       incontour=xNew<double>(nods_prime);
       for (i=0;i<nods_prime;i++) incontour[i]=1.0;
   }
 
   /*initialize thread parameters: */
   InterpFromMesh2dxThreadStruct gate;
   gate.interpolation_type = interpolation_type;
   gate.debug              = debug;
   gate.nels_data          = nels_data;
   gate.index_data         = index_data;
   gate.x_data             = x_data;
   gate.y_data             = y_data;
   gate.data               = data;
   gate.xmin               = xmin;
   gate.xmax               = xmax;
   gate.ymin               = ymin;
   gate.ymax               = ymax;
   gate.nods_prime         = nods_prime;
   gate.data_prime         = data_prime;
   gate.x_prime            = x_prime;
   gate.y_prime            = y_prime;
   gate.default_values     = default_values;
   gate.num_default_values = num_default_values;
   gate.incontour          = incontour;
 
   /*launch the thread manager with InterpFromGridToMeshxt as a core: */
   LaunchThread(InterpFromMesh2dxt,(void*)&gate,_NUMTHREADS_);
 
   /*Assign output pointers:*/
    xDelete<double>(incontour);
   *pdata_prime=data_prime;
   return 1;
}

InterpFromMesh2dxt()

void* InterpFromMesh2dxt

(

void *

vInterpFromMesh2dxThreadStruct

)

Definition at line 8 of file InterpFromMesh2dxt.cpp.

                                               {
 
   /*intermediary: */
   int     i0,i1,i,j;
   double  area,area_1,area_2,area_3;
   double  data_value;
 
   /*recover handle and gate: */
   pthread_handle                *handle      = (pthread_handle*)vpthread_handle;
   InterpFromMesh2dxThreadStruct *gate        = (InterpFromMesh2dxThreadStruct*)handle->gate;
   int                            my_thread   = handle->id;
   int                            num_threads = handle->num;
 
   /*recover parameters :*/
   int     interpolation_type      = gate->interpolation_type;
   bool    debug                   = gate->debug;
   int     nels_data               = gate->nels_data;
   int    *index_data              = gate->index_data;
   double *x_data                  = gate->x_data;
   double *y_data                  = gate->y_data;
   double *data                    = gate->data;
   double  xmin                    = gate->xmin;
   double  xmax                    = gate->xmax;
   double  ymin                    = gate->ymin;
   double  ymax                    = gate->ymax;
   int     nods_prime              = gate->nods_prime;
   IssmSeqVec<IssmPDouble>* data_prime = gate->data_prime;
   double *x_prime                 = gate->x_prime;
   double *y_prime                 = gate->y_prime;
   double *default_values          = gate->default_values;
   int     num_default_values      = gate->num_default_values;
   double *incontour               = gate->incontour;
 
   /*partition loop across threads: */
   PartitionRange(&i0,&i1,nels_data,num_threads,my_thread);
 
   /*Loop over the elements*/
   for(i=i0;i<i1;i++){
 
      /*display current iteration*/
      if (debug && my_thread==0 && fmod((double)i,(double)100)==0)
       _printf_("\r      interpolation progress: "<<setw(6)<<setprecision(2)<<double(i-i0)/double(i1-i0)*100<<"%   ");
 
      /*if there is no point inside the domain, go to next iteration*/
      if ( (x_data[index_data[3*i+0]-1]<xmin) && (x_data[index_data[3*i+1]-1]<xmin) && (x_data[index_data[3*i+2]-1]<xmin)) continue;
      if ( (x_data[index_data[3*i+0]-1]>xmax) && (x_data[index_data[3*i+1]-1]>xmax) && (x_data[index_data[3*i+2]-1]>xmax)) continue;
      if ( (y_data[index_data[3*i+0]-1]<ymin) && (y_data[index_data[3*i+1]-1]<ymin) && (y_data[index_data[3*i+2]-1]<ymin)) continue;
      if ( (y_data[index_data[3*i+0]-1]>ymax) && (y_data[index_data[3*i+1]-1]>ymax) && (y_data[index_data[3*i+2]-1]>ymax)) continue;
 
      /*get area of the current element (Jacobian = 2 * area)*/
      //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
      area=x_data[index_data[3*i+1]-1]*y_data[index_data[3*i+2]-1]-y_data[index_data[3*i+1]-1]*x_data[index_data[3*i+2]-1]
        +  x_data[index_data[3*i+0]-1]*y_data[index_data[3*i+1]-1]-y_data[index_data[3*i+0]-1]*x_data[index_data[3*i+1]-1]
        +  x_data[index_data[3*i+2]-1]*y_data[index_data[3*i+0]-1]-y_data[index_data[3*i+2]-1]*x_data[index_data[3*i+0]-1];
 
      /*loop over the prime nodes*/
      for (j=0;j<nods_prime;j++){
 
         if(incontour[j]){
 
            /*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
            area_1=((x_prime[j]-x_data[index_data[3*i+2]-1])*(y_data[index_data[3*i+1]-1]-y_data[index_data[3*i+2]-1]) 
                  -  (y_prime[j]-y_data[index_data[3*i+2]-1])*(x_data[index_data[3*i+1]-1]-x_data[index_data[3*i+2]-1]))/area;
            /*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
            area_2=((x_data[index_data[3*i+0]-1]-x_data[index_data[3*i+2]-1])*(y_prime[j]-y_data[index_data[3*i+2]-1]) 
                  - (y_data[index_data[3*i+0]-1]-y_data[index_data[3*i+2]-1])*(x_prime[j]-x_data[index_data[3*i+2]-1]))/area;
            /*Get third area coordinate = 1-area1-area2*/
            area_3=1-area_1-area_2;
 
            /*is the current point in the current element?*/
            if (area_1>=0 && area_2>=0 && area_3>=0){
 
               /*Yes ! compute the value on the point*/
               if (interpolation_type==1){
                  /*nodal interpolation*/
                  data_value=area_1*data[index_data[3*i+0]-1]+area_2*data[index_data[3*i+1]-1]+area_3*data[index_data[3*i+2]-1];
               }
               else{
                  /*element interpolation*/
                  data_value=data[i];
               }
               if (xIsNan<IssmPDouble>(data_value)){
                  if(num_default_values==1) data_value=default_values[0];
                  else data_value=default_values[j];
               }
 
               /*insert value and go to the next point*/
               data_prime->SetValue(j,data_value,INS_VAL);
            }
         }
      }
   }
   if(debug && my_thread==0)
    _printf_("\r      interpolation progress: "<<fixed<<setw(6)<<setprecision(2)<<100.<<"%  \n");
   return NULL;
}