Changeset 1189


Ignore:
Timestamp:
06/30/09 15:26:53 (16 years ago)
Author:
Mathieu Morlighem
Message:

Added x layer for InterpFromMesh*d

Location:
issm/trunk/src
Files:
6 added
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r1172 r1189  
    241241                                        ./InterpFromGridx/InterpFromGridx.cpp\
    242242                                        ./InterpFromGridx/InterpFromGridx.h\
     243                                        ./InterpFromMesh2dx/InterpFromMesh2dx.cpp\
     244                                        ./InterpFromMesh2dx/InterpFromMesh2dx.h\
     245                                        ./InterpFromMesh3dx/InterpFromMesh3dx.cpp\
     246                                        ./InterpFromMesh3dx/InterpFromMesh3dx.h\
    243247                                        ./HoleFillerx/HoleFillerx.cpp\
    244248                                        ./HoleFillerx/HoleFillerx.h\
    … …  
    521525                                        ./InterpFromGridx/InterpFromGridx.cpp\
    522526                                        ./InterpFromGridx/InterpFromGridx.h\
     527                                        ./InterpFromMesh2dx/InterpFromMesh2dx.cpp\
     528                                        ./InterpFromMesh2dx/InterpFromMesh2dx.h\
     529                                        ./InterpFromMesh3dx/InterpFromMesh3dx.cpp\
     530                                        ./InterpFromMesh3dx/InterpFromMesh3dx.h\
    523531                                        ./HoleFillerx/HoleFillerx.cpp\
    524532                                        ./HoleFillerx/HoleFillerx.h\

  • issm/trunk/src/c/issm.h

    r1172 r1189  
    2626#include "./ContourToNodesx/ContourToNodesx.h"
    2727#include "./InterpFromGridx/InterpFromGridx.h"
     28#include "./InterpFromMesh2dx/InterpFromMesh2dx.h"
     29#include "./InterpFromMesh3dx/InterpFromMesh3dx.h"
    2830#include "./HoleFillerx/HoleFillerx.h"
    2931#include "./MeshPartitionx/MeshPartitionx.h"

  • issm/trunk/src/mex/InterpFromMesh2d/InterpFromMesh2d.cpp

    r1180 r1189  
    4646
    4747        /*Intermediary*/
    48         int i,j;
    4948        int nods_data;
    5049        int nels_data;
    5150        int nods_prime;
    52         int interpolation_type;
    53         double area;
    54         double area_1,area_2,area_3;
    55         double data_value;
    56         double xmin,xmax;
    57         double ymin,ymax;
    5851
    5952        /* output: */
    … …  
    7568        FetchData((void**)&default_value,NULL,NULL,DEFAULTHANDLE,"Scalar",NULL);
    7669
    77         /* Run core computations: */
    78         //InterpFromMesh2dx( &data_mesh, index , index_rows, x, y, x_rows, data, data_rows, x_mesh, y_mesh, x_mesh_rows);
    79 
    8070        /*some checks*/
    81         if (index_data_rows<1 || x_data_rows<3 || y_data_rows<3){
    82                 throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input");
    83         }
    8471        if (x_data_rows!=y_data_rows){
    8572                throw ErrorException(__FUNCT__,"vectors x and y should have the same length!");
    … …  
    9481        nods_prime=x_prime_rows;
    9582
    96         /*figure out what kind of interpolation is needed*/
    97         if (data_rows==nods_data){
    98                 interpolation_type=1;
    99         }
    100         else if (data_rows==nels_data){
    101                 interpolation_type=2;
    102         }
    103         else{
    104                 throw ErrorException(__FUNCT__,"length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
    105         }
    106 
    107         /*Get prime mesh extrema coordinates*/
    108         xmin=x_prime[0]; xmax=x_prime[0];ymin=y_prime[0]; ymax=y_prime[0];
    109         for (i=1;i<nods_prime;i++){
    110                 if (x_prime[i]<xmin) xmin=x_prime[i];
    111                 if (x_prime[i]>xmax) xmax=x_prime[i];
    112                 if (y_prime[i]<ymin) ymin=y_prime[i];
    113                 if (y_prime[i]>ymax) ymax=y_prime[i];
    114         }
    115 
    116         /*Initialize output*/
    117         data_prime=NewVec(nods_prime);
    118         for (i=0;i<nods_prime;i++) VecSetValue(data_prime,i,default_value,INSERT_VALUES);
    119 
    120         /*Loop over the elements*/
    121         printf("\n      interpolation progress:   %5.2lf %%",0.0);
    122         for (i=0;i<nels_data;i++){
    123 
    124                 /*display current iteration*/
    125                 if (fmod(i,100)==0){
    126                         printf("\b\b\b\b\b\b\b%5.2lf %%",(double)i/nels_data*100);
    127                 }
    128 
    129                 /*if there is no point inside the domain, go to next iteration*/
    130                 if ( (x_data[(int)index_data[3*i+0]-1]<xmin) && (x_data[(int)index_data[3*i+1]-1]<xmin) && (x_data[(int)index_data[3*i+2]-1]<xmin)) continue;
    131                 if ( (x_data[(int)index_data[3*i+0]-1]>xmax) && (x_data[(int)index_data[3*i+1]-1]>xmax) && (x_data[(int)index_data[3*i+2]-1]>xmax)) continue;
    132                 if ( (y_data[(int)index_data[3*i+0]-1]<ymin) && (y_data[(int)index_data[3*i+1]-1]<ymin) && (y_data[(int)index_data[3*i+2]-1]<ymin)) continue;
    133                 if ( (y_data[(int)index_data[3*i+0]-1]>ymax) && (y_data[(int)index_data[3*i+1]-1]>ymax) && (y_data[(int)index_data[3*i+2]-1]>ymax)) continue;
    134 
    135                 /*get area of the current element (Jacobian = 2 * area)*/
    136                 //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
    137                 area=x_data[(int)index_data[3*i+1]-1]*y_data[(int)index_data[3*i+2]-1]-y_data[(int)index_data[3*i+1]-1]*x_data[(int)index_data[3*i+2]-1]+ x_data[(int)index_data[3*i+0]-1]*y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+0]-1]*x_data[(int)index_data[3*i+1]-1]+ x_data[(int)index_data[3*i+2]-1]*y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1]*x_data[(int)index_data[3*i+0]-1];
    138 
    139                 /*loop over the prime nodes*/
    140                 for (j=0;j<nods_prime;j++){
    141 
    142                         /*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
    143                         area_1=((x_prime[j]-x_data[(int)index_data[3*i+2]-1])*(y_data[(int)index_data[3*i+1]-1]-y_data[(int)index_data[3*i+2]-1]) -  (y_prime[j]-y_data[(int)index_data[3*i+2]-1])*(x_data[(int)index_data[3*i+1]-1]-x_data[(int)index_data[3*i+2]-1]))/area;
    144                         /*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
    145                         area_2=((x_data[(int)index_data[3*i+0]-1]-x_data[(int)index_data[3*i+2]-1])*(y_prime[j]-y_data[(int)index_data[3*i+2]-1]) - (y_data[(int)index_data[3*i+0]-1]-y_data[(int)index_data[3*i+2]-1])*(x_prime[j]-x_data[(int)index_data[3*i+2]-1]))/area;
    146                         /*Get third area coordinate = 1-area1-area2*/
    147                         area_3=1-area_1-area_2;
    148 
    149                         /*is the current point in the current element?*/
    150                         if (area_1>=0 && area_2>=0 && area_3>=0){
    151 
    152                                 /*Yes ! compute the value on the point*/
    153                                 if (interpolation_type==1){
    154                                         /*nodal interpolation*/
    155                                         data_value=area_1*data[(int)index_data[3*i+0]-1]+area_2*data[(int)index_data[3*i+1]-1]+area_3*data[(int)index_data[3*i+2]-1];
    156                                 }
    157                                 else{
    158                                         /*element interpolation*/
    159                                         data_value=data[i];
    160                                 }
    161                                 if isnan(data_value) data_value=default_value;
    162 
    163                                 /*insert value and go to the next point*/
    164                                 VecSetValue(data_prime,j,data_value,INSERT_VALUES);
    165                         }
    166                 }
    167         }
    168         printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
     83        /* Run core computations: */
     84        InterpFromMesh2dx(&data_prime,index_data,x_data,y_data,nods_data,nels_data,data,data_rows,x_prime,y_prime,nods_prime,default_value);
    16985
    17086        /*Write data: */

  • issm/trunk/src/mex/InterpFromMesh3d/InterpFromMesh3d.cpp

    r1180 r1189  
    5050
    5151        /*Intermediary*/
    52         int i,j;
    5352        int nods_data;
    5453        int nels_data;
    5554        int nods_prime;
    56         int interpolation_type;
    57         double area;
    58         double area_1,area_2,area_3;
    59         double zeta,bed,surface;
    60         double data_value;
    61         double xmin,xmax;
    62         double ymin,ymax;
    6355
    6456        /* output: */
    … …  
    8375
    8476        /*some checks*/
    85         if (index_data_rows<1 || x_data_rows<6 || y_data_rows<6 || z_data_rows<6){
    86                 throw ErrorException(__FUNCT__,"nothing to be done according to the mesh given in input");
    87         }
    8877        if (x_data_rows!=y_data_rows || x_data_rows!=z_data_rows){
    8978                throw ErrorException(__FUNCT__,"vectors x, y and z should have the same length!");
    … …  
    9281                throw ErrorException(__FUNCT__,"vectors x_prime, y_prime and z_prime should have the same length!");
    9382        }
    94        
    9583        /*get number of elements and number of nodes in the data*/
    9684        nels_data=index_data_rows;
    … …  
    9886        nods_prime=x_prime_rows;
    9987
    100         /*figure out what kind of interpolation is needed*/
    101         if (data_rows==nods_data){
    102                 interpolation_type=1;
    103         }
    104         else if (data_rows==nels_data){
    105                 interpolation_type=2;
    106         }
    107         else{
    108                 throw ErrorException(__FUNCT__,"length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!");
    109         }
    110 
    111         /*Get prime mesh extrema coordinates*/
    112         xmin=x_prime[0]; xmax=x_prime[0]; ymin=y_prime[0]; ymax=y_prime[0];
    113 
    114         for (i=1;i<nods_prime;i++){
    115                 if (x_prime[i]<xmin) xmin=x_prime[i];
    116                 if (x_prime[i]>xmax) xmax=x_prime[i];
    117                 if (y_prime[i]<ymin) ymin=y_prime[i];
    118                 if (y_prime[i]>ymax) ymax=y_prime[i];
    119         }
    120 
    121         /*Initialize output*/
    122         data_prime=NewVec(nods_prime);
    123         for (i=0;i<nods_prime;i++) VecSetValue(data_prime,i,default_value,INSERT_VALUES);
    124 
    125         /*Loop over the elements*/
    126         printf("\n      interpolation progress:   %5.2lf %%",0.0);
    127         for (i=0;i<nels_data;i++){
    128 
    129                 /*display current iteration*/
    130                 if (fmod(i,100)==0){
    131                         printf("\b\b\b\b\b\b\b%5.2lf %%",(double)i/nels_data*100);
    132                 }
    133 
    134                 /*if there is no point inside the domain, go to next iteration*/
    135                 if ( (x_data[(int)index_data[6*i+0]-1]<xmin) && (x_data[(int)index_data[6*i+1]-1]<xmin) && (x_data[(int)index_data[6*i+2]-1]<xmin)) continue;
    136                 if ( (x_data[(int)index_data[6*i+0]-1]>xmax) && (x_data[(int)index_data[6*i+1]-1]>xmax) && (x_data[(int)index_data[6*i+2]-1]>xmax)) continue;
    137                 if ( (y_data[(int)index_data[6*i+0]-1]<ymin) && (y_data[(int)index_data[6*i+1]-1]<ymin) && (y_data[(int)index_data[6*i+2]-1]<ymin)) continue;
    138                 if ( (y_data[(int)index_data[6*i+0]-1]>ymax) && (y_data[(int)index_data[6*i+1]-1]>ymax) && (y_data[(int)index_data[6*i+2]-1]>ymax)) continue;
    139 
    140                 /*get area of the current element (Jacobian = 2 * area)*/
    141                 //area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1;
    142                 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]+ 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]+ 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];
    143 
    144                 /*loop over the prime nodes*/
    145                 for (j=0;j<nods_prime;j++){
    146 
    147                         /*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
    148                         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]) -  (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;
    149                         /*Get second area coordinate =det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
    150                         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]) - (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;
    151                         /*Get third area coordinate = 1-area1-area2*/
    152                         area_3=1-area_1-area_2;
    153 
    154                         /*is the current point in the current 2d element?*/
    155                         if (area_1>=0 && area_2>=0 && area_3>=0){
    156 
    157                                 /*compute bottom and top height of the element at this 2d position*/
    158                                 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];
    159                                 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];
    160 
    161                                 /*Compute zeta*/
    162                                 zeta=2*(z_prime[j]-bed)/(surface-bed)-1;
    163 
    164                                 if (zeta >=-1 && zeta<=1){
    165                                         if (interpolation_type==1){
    166                                                 /*nodal interpolation*/
    167                                                 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]) + (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]);
    168                                         }
    169                                         else{
    170                                                 /*element interpolation*/
    171                                                 data_value=data[i];
    172                                         }
    173                                         if isnan(data_value) data_value=default_value;
    174 
    175                                         /*insert value and go to the next point*/
    176                                         VecSetValue(data_prime,j,data_value,INSERT_VALUES);
    177                                 }
    178                         }
    179                 }
    180         }
    181         printf("\b\b\b\b\b\b\b%5.2lf %%\n",100.0);
     88        /* Run core computations: */
     89        InterpFromMesh3dx(&data_prime,index_data,x_data,y_data,z_data,nods_data,nels_data,data,data_rows,x_prime,y_prime,z_prime,nods_prime,default_value);
    18290
    18391        /*Write data: */
Note: See TracChangeset for help on using the changeset viewer.