Changeset 5354


Ignore:
Timestamp:
08/18/10 10:05:16 (15 years ago)
Author:
Eric.Larour
Message:

New TriaSearch module.

Location:
issm/trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r5351 r5354  
    1 /*!\file InterpFromMeshToMesh2dx
     1/*!\file TriaSearchx
    22 */
    33
    4 #include "./InterpFromMeshToMesh2dx.h"
     4#include "./TriaSearchx.h"
    55
    66#include "../../shared/shared.h"
    … …  
    1313using namespace std;
    1414
    15 int InterpFromMeshToMesh2dx(double** pdata_interp,double* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
    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        
     15void TriaSearchx(double* ptria,double* index,int nel, double* x, double* y, int nods,double x0, double y0){
     16
    1817        /*Output*/
    19         double* data_interp=NULL;
     18        double tria;
    2019
    2120        /*Intermediary*/
    … …  
    2322        I2     I;
    2423        int    i,j,k;
    25         int    it;
    2624        int    i0,i1,i2;
    2725        double areacoord[3];
    … …  
    3129        int verbose=0;
    3230
    33         /*default values: */
    34         Vec    vec_incontour=NULL;
    35         double*    incontour=NULL;
    36         bool   skip_bamg=false;
    37 
    38 
    39         /*Checks*/
    40         if (data_cols<=0){
    41                 ISSMERROR("data provided has a negative number of columns");
    42         }
    43         if (data_rows!=nods_data && data_rows!=nels_data){
    44                 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);
    54         }
    55 
    56         /*Initialize output*/
    57         if (verbose) printf("Initializing output vector\n");
    58         data_interp=(double*)xmalloc(nods_interp*data_cols*sizeof(double));
    59 
    6031        // read background mesh
    6132        if (verbose) printf("Reading mesh\n");
    62         Mesh Th(index_data,x_data,y_data,nods_data,nels_data);
     33        Mesh Th(index,x,y,nods,nel);
    6334        Th.CreateSingleVertexToTriangleConnectivity();
    6435
    65         //Loop over output nodes
    66         if (verbose) printf("Loop over the nodes\n");
    67         for(i=0;i<nods_interp;i++){
    68                
    69                 /*reset skip_bamg: */
    70                 skip_bamg=false;
     36        //Get current point coordinates
     37        r.x=x0; r.y=y0;
     38        I=Th.R2ToI2(r);
    7139
    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                         }
    78                 }
     40        //Find triangle holding r/I
     41        Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
    7942
    80                 if(skip_bamg==false){
    81 
    82                         //Get current point coordinates
    83                         r.x=x_interp[i]; r.y=y_interp[i];
    84                         I2 I=Th.R2ToI2(r);
    85 
    86                         //Find triangle holding r/I
    87                         Triangle &tb=*Th.TriangleFindFromCoord(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.GetId(tb[0]);
    97                                 i1=Th.GetId(tb[1]);
    98                                 i2=Th.GetId(tb[2]);
    99                                 //triangle number
    100                                 it=Th.GetId(tb);
    101                         }
    102                         //external point
    103                         else {
    104                                 //Get closest adjacent triangle (inside the mesh)
    105                                 AdjacentTriangle 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.GetId(tc[0]);
    114                                 i1=Th.GetId(tc[1]);
    115                                 i2=Th.GetId(tc[2]);
    116                                 //triangle number
    117                                 it=Th.GetId(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                                 }
    132                         }
    133                 }
    134                 else{
    135                         if(num_default_values==1) data_interp[i]=default_values[0];
    136                         else data_interp[i]=default_values[i];
    137                 }
     43        // internal point
     44        if (tb.det>0){
     45                //Area coordinate
     46                areacoord[0]= (double) dete[0]/ tb.det;
     47                areacoord[1]= (double) dete[1] / tb.det;
     48                areacoord[2]= (double) dete[2] / tb.det;
     49                //3 vertices of the triangle
     50                i0=Th.GetId(tb[0]);
     51                i1=Th.GetId(tb[1]);
     52                i2=Th.GetId(tb[2]);
     53                //triangle number
     54                tria=(double)Th.GetId(tb);
    13855        }
    139 
     56        //external point
     57        else tria=NAN;
     58       
    14059        /*Assign output pointers:*/
    141         if (verbose) printf("Assigning output\n");
    142         *pdata_interp=data_interp;
    143 
    144         /*No error return*/
    145         return 1;
     60        *ptria=tria;
    14661}

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

    r5351 r5354  
    1 /*!\file:  InterpFromMeshToMesh2dx.h
     1/*!\file:  TriaSearchx.h
    22 * \brief header file for Bamg module
    33 */
    44
    5 #ifndef _INTERPFROMMESHTOMESH2DX_H
    6 #define _INTERPFROMMESHTOMESH2DX_H
     5#ifndef _TRIASEARCHX_H
     6#define _TRIASEARCHX_H
    77
    88#include "../../objects/objects.h"
    99
    1010/* local prototypes: */
    11 int 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_values,int num_default_values,Contour** contours, int numcontours );
     11void TriaSearchx(double* ptria,double* index,int nel, double* x, double* y, int nods,double x0, double y0);
    1312
    1413#endif

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

    r5249 r5354  
    6060                                SystemMatrices\
    6161                                Test\
     62                                TriaSearch\
    6263                                TriMesh\
    6364                                TriMeshNoDensity\
    … …  
    260261                          SystemMatrices/SystemMatrices.h
    261262
     263TriaSearch_SOURCES = TriaSearch/TriaSearch.cpp\
     264                          TriaSearch/TriaSearch.h
     265
    262266TriMesh_SOURCES = TriMesh/TriMesh.cpp\
    263267                          TriMesh/TriMesh.h

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

    r5348 r5354  
    1 /*\file InterpFromMeshToMesh2d.c
    2  *\brief: bamg module.
     1/*\file TriaSearch.c
     2 *\brief: TriaSearch module. See TriaSearchx for more details.
    33 */
    4 #include "./InterpFromMeshToMesh2d.h"
     4#include "./TriaSearch.h"
    55
    66void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){
    … …  
    1010        /*input: */
    1111        double* index=NULL;
    12         int     index_cols;
     12        int     nel;
     13        int     dummy;
    1314
    14         double* x_data=NULL;
    15         int     x_data_rows;
     15        double* x=NULL;
     16        int     nods;
    1617
    17         double* y_data=NULL;
    18         int     y_data_rows;
     18        double* y=NULL;
    1919
    20         double* data=NULL;
    21         int     data_rows;
    22         int     data_cols;
    23 
    24         double* x_interp=NULL;
    25         double* y_interp=NULL;
    26 
    27         int     x_interp_rows;
    28         int     y_interp_rows;
    29 
    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;
    38 
    39         /*Intermediary*/
    40         int nods_data;
    41         int nels_data;
    42         int nods_interp;
    43         int verbose=0;
    44 
     20        double  x0,y0;
     21       
    4522        /* output: */
    46         double* data_interp=NULL;
     23        double  tria;
    4724
    4825        /*Boot module: */
    … …  
    5027
    5128        /*checks on arguments on the matlab side: */
    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         }
     29        CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&TriaSearchUsage);
    6030
    6131        /*Input datasets: */
    62         if (verbose) printf("Fetching inputs\n");
    63         FetchData(&index,&nels_data,&index_cols,INDEXHANDLE);
    64         FetchData(&x_data,&x_data_rows,NULL,XHANDLE);
    65         FetchData(&y_data,&y_data_rows,NULL,YHANDLE);
    66         FetchData(&data,&data_rows,&data_cols,DATAHANDLE);
    67         FetchData(&x_interp,&x_interp_rows,NULL,XINTERPHANDLE);
    68         FetchData(&y_interp,&y_interp_rows,NULL,YINTERPHANDLE);
     32        FetchData(&index,&nel,&dummy,INDEXHANDLE);
     33        FetchData(&x,&nods,XHANDLE);
     34        FetchData(&y,&nods,YHANDLE);
     35        FetchData(&x0,X0HANDLE);
     36        FetchData(&y0,Y0HANDLE);
    6937
    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         }
    109 
    110         /*some checks*/
    111         if (verbose) printf("Checking inputs\n");
    112         if (x_data_rows!=y_data_rows){
    113                 ISSMERROR("vectors x and y should have the same length!");
    114         }
    115         if (x_interp_rows!=y_interp_rows){
    116                 ISSMERROR("vectors x_interp and y_interp should have the same length!");
    117         }
    118         if (index_cols!=3){
    119                 ISSMERROR("index should have 3 columns (input provided has %i columns)",index_cols);
    120         }
    121 
    122         /*get number of elements and number of nodes in the data*/
    123         nods_data=x_data_rows;
    124         nods_interp=x_interp_rows;
     38        /* Echo: {{{1*/
     39        //printf("(x0,y0)=(%g,%g)\n",x0,y0);
     40        /*}}}*/
    12541
    12642        /* Run core computations: */
    127         if (verbose) printf("Call core\n");
    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);
     43        TriaSearchx(&tria,index,nel,x,y,nods,x0,y0);
     44
     45        /* c to matlab: */
     46        tria++;
    12947
    13048        /*Write data: */
    131         WriteData(DATAINTERP,data_interp,nods_interp,data_cols);
     49        WriteData(TRIA,tria);
    13250
    13351        /*end module: */
    … …  
    13553}
    13654
    137 void InterpFromMeshToMesh2dUsage(void)
     55void TriaSearchUsage(void)
    13856{
    139         _printf_("INTERFROMMESHTOMESH2D - interpolation from a 2d triangular mesh onto a list of point\n");
    140         _printf_("\n");
    141         _printf_("   This function is a multi-threaded mex file that interpolates a field\n");
    142         _printf_("   defined on a triangular mesh onto a list of point\n");
     57        _printf_("TriaSearch- find triangle holding a point (x0,y0) in a mesh\n");
    14358        _printf_("\n");
    14459        _printf_("   Usage:\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");
    147         _printf_("\n");
    148         _printf_("      index: index of the mesh where data is defined\n");
    149         _printf_("      x,y: coordinates of the nodes where data is defined\n");
    150         _printf_("      data: matrix holding the data to be interpolated onto the mesh. (one column per field)\n");
    151         _printf_("      x_interp,y_interp: coordinates of the points onto which we interpolate.\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");
    155         _printf_("      data_interp: vector of mesh interpolated data.\n");
    156         _printf_("\n");
    157         _printf_("   Example:\n");
    158         _printf_("      load('temperature.mat');\n");
    159         _printf_("      md.temperature=InterpFromMeshToMesh2d(index,x,y,temperature,md.x,md.y);\n");
    160         _printf_("      md.temperature=InterpFromMeshToMesh2d(index,x,y,temperature,md.x,md.y,253,'Contour.exp');\n");
     60        _printf_("         tria=TriaSearch(index,x,y,x0,y0);\n");
     61        _printf_("      index,x,y: mesh triangulatrion\n");
     62        _printf_("      x0,y0: coordinates of the point for which we are trying to find a triangle\n");
    16163        _printf_("\n");
    16264}

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

    r5348 r5354  
    1 /*!\file InterpFromMeshToMesh2d.h
    2  * \brief: prototype for Data Interpolation mex module.
     1/*!\file TriaSearch.h
    32 */
    43
    5 #ifndef _INTERPFROMMESHTOMESH2d_H
    6 #define _INTERPFROMMESHTOMESH2d_H
     4#ifndef _TRIASEARCH_H
     5#define _TRIASEARCH_H
    76
    87/* local prototypes: */
    9 void InterpFromMeshToMesh2dUsage(void);
     8void TriaSearchUsage(void);
    109
    1110#include "../../c/modules/modules.h"
    … …  
    1413
    1514#undef __FUNCT__
    16 #define __FUNCT__  "InterpFromMeshToMesh2d"
     15#define __FUNCT__  "TriaSearch"
    1716
    1817
    … …  
    2120#define XHANDLE prhs[1]
    2221#define YHANDLE prhs[2]
    23 #define DATAHANDLE prhs[3]
    24 #define XINTERPHANDLE prhs[4]
    25 #define YINTERPHANDLE prhs[5]
    26 #define DEFAULTHANDLE prhs[6]
    27 #define FILENAME prhs[7]
     22#define X0HANDLE prhs[3]
     23#define Y0HANDLE prhs[4]
    2824
    2925/* serial output macros: */
    30 #define DATAINTERP (mxArray**)&plhs[0]
     26#define TRIA (mxArray**)&plhs[0]
    3127
    3228/* serial arg counts: */
    … …  
    3430#define NLHS  1
    3531
     32#undef NRHS
     33#define NRHS  5
     34
    3635#endif
Note: See TracChangeset for help on using the changeset viewer.