Changeset 5032


Ignore:
Timestamp:
08/06/10 14:48:41 (15 years ago)
Author:
Eric.Larour
Message:

Added InterpFromMesh2d module, from ISSM2.0, to go around some of Bamg
error messages. We'll get rid of it later, when bugs are addressed.

Location:
issm/trunk/src
Files:
7 added
7 edited

Legend:

Unmodified
Added
Removed

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

    r4983 r5032  
    479479                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
    480480                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
     481                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\
     482                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\
    481483                                        ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\
    482484                                        ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h\
    … …  
    10151017                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.cpp\
    10161018                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshx.h\
     1019                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.cpp\
     1020                                        ./modules/InterpFromMesh2dx/InterpFromMesh2dx.h\
    10171021                                        ./modules/InterpFromGridToMeshx/InterpFromGridToMeshxt.cpp\
    10181022                                        ./modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp\

  • issm/trunk/src/c/modules/InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.cpp

    r5016 r5032  
    88#include "../../toolkits/toolkits.h"
    99#include "../../objects/objects.h"
     10#include "../modules.h"
    1011
    1112using namespace bamg;
    … …  
    1314
    1415int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
    15                         double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value){
    16 
     16                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double* default_values,int num_default_values, Contour** contours, int numcontours){
     17       
    1718        /*Output*/
    1819        double* data_interp=NULL;
    … …  
    3031        int verbose=0;
    3132
     33        /*default values: */
     34        Vec    vec_incontour=NULL;
     35        double*    incontour=NULL;
     36        bool   skip_bamg=false;
     37
     38
    3239        /*Checks*/
    3340        if (data_cols<=0){
    … …  
    3643        if (data_rows!=nods_data && data_rows!=nels_data){
    3744                ISSMERROR("data provided should have either %i or %i lines (not %i)",nods_data,nels_data,data_rows);
     45        }
     46        if((num_default_values) && (data_cols>1)){
     47                ISSMERROR("data provided can only have 1 column if a default value is provided");
     48        }
     49       
     50        /*If default values supplied, figure out which nodes are inside the contour, including the border of the contour: */
     51        if(num_default_values){
     52                ContourToNodesx( &vec_incontour,x_interp,y_interp,nods_interp,contours,numcontours,1);
     53                VecToMPISerial(&incontour,vec_incontour);
    3854        }
    3955
    … …  
    5066        if (verbose) printf("Loop over the nodes\n");
    5167        for(i=0;i<nods_interp;i++){
     68               
     69                /*reset skip_bamg: */
     70                skip_bamg=false;
    5271
    53                 //Get current point coordinates
    54                 r.x=x_interp[i]; r.y=y_interp[i];
    55                 I2 I=Th.toI2(r);
    56 
    57                 //Find triangle holding r/I
    58                 Triangle &tb=*Th.FindTriangleContaining(I,dete);
    59 
    60                 // internal point
    61                 if (tb.det>0){
    62                         //Area coordinate
    63                         areacoord[0]= (double) dete[0]/ tb.det;
    64                         areacoord[1]= (double) dete[1] / tb.det;
    65                         areacoord[2]= (double) dete[2] / tb.det;
    66                         //3 vertices of the triangle
    67                         i0=Th.Number(tb[0]);
    68                         i1=Th.Number(tb[1]);
    69                         i2=Th.Number(tb[2]);
    70                         //triangle number
    71                         it=Th.Number(tb);
     72                /*figure out if we should skip bamg logic: */
     73                if(num_default_values){
     74                        if(!incontour[i]){
     75                                /*This node is not inside the contour. Skip Bamg logic and apply default value.: */
     76                                skip_bamg=true;
     77                        }
    7278                }
    7379
    74                 //external point
    75                 else {
    76                         //Get closest adjacent triangle (inside the mesh)
    77                         TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
    78                         int k=ta;
    79                         Triangle &tc=*(Triangle*)ta;
    80                         //Area coordinate
    81                         areacoord[VerticesOfTriangularEdge[k][1]] = aa;
    82                         areacoord[VerticesOfTriangularEdge[k][0]] = bb;
    83                         areacoord[OppositeVertex[k]] = 1 - aa -bb;
    84                         //3 vertices of the triangle
    85                         i0=Th.Number(tc[0]);
    86                         i1=Th.Number(tc[1]);
    87                         i2=Th.Number(tc[2]);
    88                         //triangle number
    89                         it=Th.Number(tc);
    90                 }
     80                if(skip_bamg==false){
    9181
    92                 /*Last step, P1 interpolation*/
    93                 if (data_rows==nods_data){
    94                         for (j=0;j<data_cols;j++){
    95                                 data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j];
     82                        //Get current point coordinates
     83                        r.x=x_interp[i]; r.y=y_interp[i];
     84                        I2 I=Th.toI2(r);
     85
     86                        //Find triangle holding r/I
     87                        Triangle &tb=*Th.FindTriangleContaining(I,dete);
     88
     89                        // internal point
     90                        if (tb.det>0){
     91                                //Area coordinate
     92                                areacoord[0]= (double) dete[0]/ tb.det;
     93                                areacoord[1]= (double) dete[1] / tb.det;
     94                                areacoord[2]= (double) dete[2] / tb.det;
     95                                //3 vertices of the triangle
     96                                i0=Th.Number(tb[0]);
     97                                i1=Th.Number(tb[1]);
     98                                i2=Th.Number(tb[2]);
     99                                //triangle number
     100                                it=Th.Number(tb);
     101                        }
     102                        //external point
     103                        else {
     104                                //Get closest adjacent triangle (inside the mesh)
     105                                TriangleAdjacent ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
     106                                int k=ta;
     107                                Triangle &tc=*(Triangle*)ta;
     108                                //Area coordinate
     109                                areacoord[VerticesOfTriangularEdge[k][1]] = aa;
     110                                areacoord[VerticesOfTriangularEdge[k][0]] = bb;
     111                                areacoord[OppositeVertex[k]] = 1 - aa -bb;
     112                                //3 vertices of the triangle
     113                                i0=Th.Number(tc[0]);
     114                                i1=Th.Number(tc[1]);
     115                                i2=Th.Number(tc[2]);
     116                                //triangle number
     117                                it=Th.Number(tc);
     118                        }
     119                       
     120                        if (data_rows==nods_data){
     121                                for (j=0;j<data_cols;j++){
     122                                        data_interp[i*data_cols+j]=areacoord[0]*data[data_cols*i0+j]+areacoord[1]*data[data_cols*i1+j]+areacoord[2]*data[data_cols*i2+j];
     123                                }
     124                        }
     125                        else{
     126                                for (j=0;j<data_cols;j++){
     127                                        if (it<0 || it>=nels_data){
     128                                                ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data);
     129                                        }
     130                                        data_interp[i*data_cols+j]=data[data_cols*it+j];
     131                                }
    96132                        }
    97133                }
    98134                else{
    99                         for (j=0;j<data_cols;j++){
    100                                 if (it<0 || it>=nels_data){
    101                                         ISSMERROR("Triangle number %i not in [0 %i], because not correctly implemented yet... interpolate on grid first",it,nels_data);
    102                                 }
    103                                 data_interp[i*data_cols+j]=data[data_cols*it+j];
    104                         }
     135                        if(num_default_values==1) data_interp[i]=default_values[0];
     136                        else data_interp[i]=default_values[i];
    105137                }
    106 
    107138        }
    108139

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

    r3913 r5032  
    1010/* local prototypes: */
    1111int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
    12                         double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double default_value);
     12                        double* data,int data_rows,int data_cols,double* x_interp,double* y_interp,int nods_interp,double* default_values,int num_default_values,Contour** contours, int numcontours );
    1313
    1414#endif

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

    r4983 r5032  
    3030#include "./InputDuplicatex/InputDuplicatex.h"
    3131#include "./InputScalex/InputScalex.h"
     32#include "./InterpFromMesh2dx/InterpFromMesh2dx.h"
    3233#include "./InterpFromGridToMeshx/InterpFromGridToMeshx.h"
    3334#include "./InterpFromMeshToMesh2dx/InterpFromMeshToMesh2dx.h"

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

    r4981 r5032  
    55
    66void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
     7
     8        int i;
    79
    810        /*input: */
    … …  
    2628        int     y_interp_rows;
    2729
    28         double default_value;
     30        double* default_values=NULL;
     31        int     num_default_values=0;
     32
     33        //contours
     34        mxArray*  matlabstructure=NULL;
     35        Contour** contours=NULL;
     36        int       numcontours;
     37        Contour*  contouri=NULL;
    2938
    3039        /*Intermediary*/
    … …  
    4150
    4251        /*checks on arguments on the matlab side: */
    43         CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&InterpFromMeshToMesh2dUsage);
     52        if(nlhs!=NLHS){
     53                InterpFromMeshToMesh2dUsage();
     54                ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
     55        }
     56        if((nrhs!=6) & (nrhs!=8)){
     57                InterpFromMeshToMesh2dUsage();
     58                ISSMERROR("InterpFromMeshToMesh2dUsage usage error");
     59        }
    4460
    4561        /*Input datasets: */
    … …  
    5167        FetchData(&x_interp,&x_interp_rows,NULL,XINTERPHANDLE);
    5268        FetchData(&y_interp,&y_interp_rows,NULL,YINTERPHANDLE);
    53         FetchData(&default_value,DEFAULTHANDLE);
     69
     70        if(nrhs==8){
     71               
     72                /*Call expread on filename to build a contour array in the matlab workspace: */
     73                mexCallMATLAB( 1, &matlabstructure, 1, (mxArray**)&FILENAME, "expread");
     74
     75                /*default values: */
     76                FetchData(&default_values,&num_default_values,DEFAULTHANDLE);
     77
     78                /*contours: */
     79                numcontours=mxGetNumberOfElements(matlabstructure);
     80                contours=(Contour**)xmalloc(numcontours*sizeof(Contour*));
     81                for(i=0;i<numcontours;i++){
     82                        //allocate
     83                        contouri=(Contour*)xmalloc(sizeof(Contour));
     84                        //retrieve dimension of this contour.
     85                        contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));
     86                        //set pointers.
     87                        contouri->x=mxGetPr(mxGetField(matlabstructure,i,"x"));
     88                        contouri->y=mxGetPr(mxGetField(matlabstructure,i,"y"));
     89                        *(contours+i)=contouri;
     90                }
     91
     92                /* Debugging of contours :{{{1*/
     93                /*for(i=0;i<numcontours;i++){
     94                  printf("\nContour echo: contour number  %i / %i\n",i+1,numcontours);
     95                  contouri=*(contours+i);
     96                  printf("   Number of grids %i\n",contouri->nods);
     97                  for (j=0;j<contouri->nods;j++){
     98                  printf("   %lf %lf\n",*(contouri->x+j),*(contouri->y+j));
     99                  }
     100                  }*/
     101                /*}}}*/
     102        }
     103        else{
     104                num_default_values=0;
     105                default_values=NULL;
     106                numcontours=0;
     107                contours=NULL;
     108        }
    54109
    55110        /*some checks*/
    … …  
    71126        /* Run core computations: */
    72127        if (verbose) printf("Call core\n");
    73         InterpFromMeshToMesh2dx(&data_interp,index,x_data,y_data,nods_data,nels_data,data,data_rows,data_cols,x_interp,y_interp,nods_interp,default_value);
     128        InterpFromMeshToMesh2dx(&data_interp,index,x_data,y_data,nods_data,nels_data,data,data_rows,data_cols,x_interp,y_interp,nods_interp,default_values,num_default_values,contours,numcontours);
    74129
    75130        /*Write data: */
    … …  
    88143        _printf_("\n");
    89144        _printf_("   Usage:\n");
    90         _printf_("      data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value);\n");
     145        _printf_("         data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp);\n");
     146        _printf_("      or data_interp=InterpFromMeshToMesh2d(index,x,y,data,x_interp,y_interp,default_value,contourname);\n");
    91147        _printf_("\n");
    92148        _printf_("      index: index of the mesh where data is defined\n");
    … …  
    94150        _printf_("      data: matrix holding the data to be interpolated onto the mesh. (one column per field)\n");
    95151        _printf_("      x_interp,y_interp: coordinates of the points onto which we interpolate.\n");
    96         _printf_("      default_value: default value if no data is found (holes).\n");
     152        _printf_("      if default_value and contourname not specified: linear interpolation will happen on all x_interp,y_interp.\n");
     153        _printf_("      if (default_value,contourname) specified: linear interpolation will happen on all x_interp,y_interp inside the contour, default value will be adopted on the rest of the mesh.\n");
     154        _printf_("      note that default_value is either a scalar, or a vector of size  length(x_interp)\n");
    97155        _printf_("      data_interp: vector of mesh interpolated data.\n");
    98156        _printf_("\n");

  • issm/trunk/src/mex/InterpFromMeshToMesh2d/InterpFromMeshToMesh2d.h

    r4236 r5032  
    2525#define YINTERPHANDLE prhs[5]
    2626#define DEFAULTHANDLE prhs[6]
     27#define FILENAME prhs[7]
    2728
    2829/* serial output macros: */
    … …  
    3233#undef NLHS
    3334#define NLHS  1
    34 #undef NRHS
    35 #define NRHS  7
    3635
    3736#endif

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

    r5004 r5032  
    3535                                InterpFromMeshToMesh3d \
    3636                                InterpFromMeshToGrid \
     37                                InterpFromMesh2d \
    3738                                MassFlux\
    3839                                Mergesolutionfromftog\
    … …  
    191192                                                                        InterpFromMeshToGrid/InterpFromMeshToGrid.h
    192193
     194InterpFromMesh2d_SOURCES = InterpFromMesh2d/InterpFromMesh2d.cpp\
     195                                                                        InterpFromMesh2d/InterpFromMesh2d.h
     196
    193197AverageFilter_SOURCES = AverageFilter/AverageFilter.cpp\
    194198                          AverageFilter/AverageFilter.h
Note: See TracChangeset for help on using the changeset viewer.