Changeset 12600


Ignore:
Timestamp:
07/03/12 08:39:17 (13 years ago)
Author:
Mathieu Morlighem
Message:

Simplified InterpFromMeshToMesh2dx by adding options, and index is now always int

Location:
issm/trunk-jpl/src/c
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/modules/BamgConvertMeshx/BamgConvertMeshx.cpp

    r12511 r12600  
    1212using namespace std;
    1313
    14 int BamgConvertMeshx(BamgMesh* bamgmesh,BamgGeom* bamggeom,double* index,double* x,double* y,int nods,int nels){
    15 
    16         /*Intermediary*/
    17         int i,j,k;
    18         int verbose=0;
    19         int noerr=1;
     14int BamgConvertMeshx(BamgMesh* bamgmesh,BamgGeom* bamggeom,int* index,double* x,double* y,int nods,int nels){
    2015
    2116        /*Options*/
    22         BamgOpts* bamgopts=NULL;
    23         bamgopts=new BamgOpts();
     17        BamgOpts* bamgopts=new BamgOpts();
    2418
    25         // read mesh
    26         if(verbose) _printLine_("Reading mesh");
     19        /*read mesh*/
    2720        Mesh Th(index,x,y,nods,nels);
    2821
    29         //write mesh and geometry
    30         if (verbose) _printLine_("Write Geometry");
     22        /*write mesh and geometry*/
    3123        Th.Gh.WriteGeometry(bamggeom,bamgopts);
    32         if (verbose) _printLine_("Write Mesh");
    3324        Th.WriteMesh(bamgmesh,bamgopts);
    3425
    35         //clean up
     26        /*clean up and return*/
    3627        delete bamgopts;
    37 
    38         /*No error return*/
    39         return noerr;
     28        return 1;
    4029
    4130}

  • issm/trunk-jpl/src/c/modules/BamgConvertMeshx/BamgConvertMeshx.h

    r3913 r12600  
    99
    1010/* local prototypes: */
    11 int BamgConvertMeshx(BamgMesh* bamgmesh,BamgGeom* bamggeom,double* index,double* x,double* y,int nods,int nels);
     11int BamgConvertMeshx(BamgMesh* bamgmesh,BamgGeom* bamggeom,int* index,double* x,double* y,int nods,int nels);
    1212
    1313#endif

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

    r12511 r12600  
    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, DataSet* contours){
     15int InterpFromMeshToMesh2dx(double** pdata_interp,int* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
     16                        double* data,int M_data,int N_data,double* x_interp,double* y_interp,int N_interp,Options* options){
    1717       
    1818        /*Output*/
    … …  
    2020
    2121        /*Intermediary*/
     22        bool   isdefault;
     23        double defaultvalue;
    2224        R2     r;
    2325        I2     I;
    … …  
    2931        double data_value;
    3032        Icoor2 dete[3];
    31         int verbose=0;
    32 
    33         /*default values: */
    34         Vector*    vec_incontour=NULL;
    35         double*    incontour=NULL;
    36         bool   skip_bamg=false;
    3733
    3834        /*Checks*/
    39         if (data_cols<=0){
    40                 _error2_("data provided has a negative number of columns");
     35        if (M_data!=nods_data && M_data!=nels_data){
     36                _error2_("data provided should have either " << nods_data << " or " << nels_data << " lines (not " << M_data << ")");
    4137        }
    42         if (data_rows!=nods_data && data_rows!=nels_data){
    43                 _error2_("data provided should have either " << nods_data << " or " << nels_data << " lines (not " << data_rows << ")");
     38
     39        /*Get default*/
     40        if(options->GetOption("default")){
     41                isdefault=true;
     42                options->Get(&defaultvalue,"default");
    4443        }
    45         if((num_default_values) && (data_cols>1)){
    46                 _error2_("data provided can only have 1 column if a default value is provided");
    47         }
    48        
    49         /*If default values supplied, figure out which nodes are inside the contour, including the border of the contour: */
    50         if(num_default_values){
    51                 ContourToNodesx( &vec_incontour,x_interp,y_interp,nods_interp,contours,1);
    52                 incontour=vec_incontour->ToMPISerial();
     44        else{
     45                isdefault=false;
    5346        }
    5447
    5548        /*Initialize output*/
    56         if (verbose) _printLine_("Initializing output vector");
    57         data_interp=xNew<double>(nods_interp*data_cols);
     49        data_interp=xNew<double>(N_interp*N_data);
    5850
    59         // read background mesh
    60         if (verbose) _printLine_("Reading mesh");
     51        /*read background mesh*/
    6152        Mesh Th(index_data,x_data,y_data,nods_data,nels_data);
    6253        Th.CreateSingleVertexToTriangleConnectivity();
    6354
    64         //Loop over output nodes
    65         if (verbose) _printLine_("Loop over the nodes");
    66         for(i=0;i<nods_interp;i++){
    67                
    68                 /*reset skip_bamg: */
    69                 skip_bamg=false;
     55        /*Loop over output nodes*/
     56        for(i=0;i<N_interp;i++){
    7057
    71                 /*figure out if we should skip bamg logic: */
    72                 if(num_default_values){
    73                         if(!incontour[i]){
    74                                 /*This node is not inside the contour. Skip Bamg logic and apply default value.: */
    75                                 skip_bamg=true;
     58                //Get current point coordinates
     59                r.x=x_interp[i]; r.y=y_interp[i];
     60                I2 I=Th.R2ToI2(r);
     61
     62                //Find triangle holding r/I
     63                Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
     64
     65                // internal point
     66                if (tb.det>0){
     67                        //Area coordinate
     68                        areacoord[0]= (double) dete[0]/tb.det;
     69                        areacoord[1]= (double) dete[1]/tb.det;
     70                        areacoord[2]= (double) dete[2]/tb.det;
     71                        //3 vertices of the triangle
     72                        i0=Th.GetId(tb[0]);
     73                        i1=Th.GetId(tb[1]);
     74                        i2=Th.GetId(tb[2]);
     75                        //triangle number
     76                        it=Th.GetId(tb);
     77                }
     78                //external point
     79                else{
     80                        if(isdefault){
     81                                for(j=0;j<N_data;j++) data_interp[i*N_data+j]=defaultvalue;
     82                                continue;
    7683                        }
    77                 }
    78 
    79                 if(skip_bamg==false){
    80 
    81                         //Get current point coordinates
    82                         r.x=x_interp[i]; r.y=y_interp[i];
    83                         I2 I=Th.R2ToI2(r);
    84 
    85                         //Find triangle holding r/I
    86                         Triangle &tb=*Th.TriangleFindFromCoord(I,dete);
    87 
    88                         // internal point
    89                         if (tb.det>0){
    90                                 //Area coordinate
    91                                 areacoord[0]= (double) dete[0]/tb.det;
    92                                 areacoord[1]= (double) dete[1]/tb.det;
    93                                 areacoord[2]= (double) dete[2]/tb.det;
    94                                 //3 vertices of the triangle
    95                                 i0=Th.GetId(tb[0]);
    96                                 i1=Th.GetId(tb[1]);
    97                                 i2=Th.GetId(tb[2]);
    98                                 //triangle number
    99                                 it=Th.GetId(tb);
    100                         }
    101                         //external point
    102                         else {
     84                        else{
    10385                                //Get closest adjacent triangle (inside the mesh)
    10486                                AdjacentTriangle ta=CloseBoundaryEdge(I,&tb,aa,bb).Adj();
    … …  
    11698                                it=Th.GetId(tc);
    11799                        }
    118                        
    119                         if (data_rows==nods_data){
    120                                 for (j=0;j<data_cols;j++){
    121                                         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];
    122                                 }
    123                         }
    124                         else{
    125                                 for (j=0;j<data_cols;j++){
    126                                         if (it<0 || it>=nels_data){
    127                                                 _error2_("Triangle number " << it << " not in [0 " << nels_data << "], because not correctly implemented yet... interpolate on grid first");
    128                                         }
    129                                         data_interp[i*data_cols+j]=data[data_cols*it+j];
    130                                 }
     100                }
     101
     102                if (M_data==nods_data){
     103                        for (j=0;j<N_data;j++){
     104                                data_interp[i*N_data+j]=areacoord[0]*data[N_data*i0+j]+areacoord[1]*data[N_data*i1+j]+areacoord[2]*data[N_data*i2+j];
    131105                        }
    132106                }
    133107                else{
    134                         if(num_default_values==1) data_interp[i]=default_values[0];
    135                         else data_interp[i]=default_values[i];
     108                        for (j=0;j<N_data;j++){
     109                                if (it<0 || it>=nels_data){
     110                                        _error2_("Triangle number " << it << " not in [0 " << nels_data << "], report bug to developers");
     111                                }
     112                                data_interp[i*N_data+j]=data[N_data*it+j];
     113                        }
    136114                }
    137115        }
    138116
    139117        /*Assign output pointers:*/
    140         if (verbose) _printLine_("Assigning output");
    141118        *pdata_interp=data_interp;
    142119

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

    r12121 r12600  
    66#define _INTERPFROMMESHTOMESH2DX_H
    77
    8 #include "../../objects/objects.h"
     8class Options;
    99
    10 /* 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,DataSet* contours);
     10int InterpFromMeshToMesh2dx(double** pdata_interp,int* index_data,double* x_data,double* y_data,int nods_data,int nels_data,
     11                        double* data,int M_data,int N_data,double* x_interp,double* y_interp,int N_interp,Options* options);
    1312
    1413#endif

  • issm/trunk-jpl/src/c/modules/Krigingx/Krigingx.cpp

    r12595 r12600  
    6262                double *data  = NULL;
    6363                int    *index = NULL;
    64                 double *indexd= NULL;
    6564
    6665                observations->ObservationList(&x,&y,&data,&nobs);
    … …  
    6867                _printLine_("Generation Delaunay Triangulation");
    6968                BamgTriangulatex(&index,&nel,x,y,nobs);
    70                 indexd =xNewZeroInit<double>(nel*3);
    71                 for(int i=0;i<nel*3;i++) indexd[i]=(double)index[i];
    72                 xDelete<int>(index);
    7369
    7470                _printLine_("Interpolating");
    7571                xDelete<double>(predictions);
    76                 InterpFromMeshToMesh2dx(&predictions,indexd,x,y,nobs,nel,data,nobs,1,x_interp,y_interp,n_interp,NULL,0,new DataSet());
     72                InterpFromMeshToMesh2dx(&predictions,index,x,y,nobs,nel,data,nobs,1,x_interp,y_interp,n_interp,options);
    7773                xDelete<double>(x);
    7874                xDelete<double>(y);

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

    r12450 r12600  
    1313using namespace std;
    1414
    15 void TriaSearchx(double** ptria,double* index,int nel, double* x, double* y, int nods,double* x0, double* y0,int numberofnodes){
     15void TriaSearchx(double** ptria,int* index,int nel, double* x, double* y, int nods,double* x0, double* y0,int numberofnodes){
    1616
    1717        /*Output*/

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

    r8303 r12600  
    99
    1010/* local prototypes: */
    11 void TriaSearchx(double** ptria,double* index,int nel, double* x, double* y, int nods,double* x0, double* y0,int numberofnodes);
     11void TriaSearchx(double** ptria,int* index,int nel, double* x, double* y, int nods,double* x0, double* y0,int numberofnodes);
    1212
    1313#endif

  • issm/trunk-jpl/src/c/objects/Bamg/Mesh.cpp

    r12574 r12600  
    4141        }
    4242        /*}}}*/
    43         /*FUNCTION Mesh::Mesh(double* index,double* x,double* y,int nods,int nels){{{*/
    44         Mesh::Mesh(double* index,double* x,double* y,int nods,int nels):Gh(*(new Geometry())),BTh(*this){
     43        /*FUNCTION Mesh::Mesh(int* index,double* x,double* y,int nods,int nels){{{*/
     44        Mesh::Mesh(int* index,double* x,double* y,int nods,int nels):Gh(*(new Geometry())),BTh(*this){
    4545
    4646                Init(0);
    … …  
    265265
    266266        /*IO*/
    267         /*FUNCTION Mesh::ReadMesh(double* index,double* x,double* y,int nods,int nels){{{*/
    268         void Mesh::ReadMesh(double* index,double* x,double* y,int nods,int nels){
     267        /*FUNCTION Mesh::ReadMesh(int* index,double* x,double* y,int nods,int nels){{{*/
     268        void Mesh::ReadMesh(int* index,double* x,double* y,int nods,int nels){
    269269
    270270                double Hmin = HUGE_VAL;// the infinie value

  • issm/trunk-jpl/src/c/objects/Bamg/Mesh.h

    r12218 r12600  
    5757                        //Constructors/Destructors
    5858                        Mesh(BamgGeom* bamggeom,BamgMesh* bamgmesh,BamgOpts* bamgopts);
    59                         Mesh(double* index,double* x,double* y,int nods,int nels);/*MeshConvert*/
     59                        Mesh(int* index,double* x,double* y,int nods,int nels);/*MeshConvert*/
    6060                        Mesh(double* x,double* y,int nods); /*BamgTriangulate*/
    6161                        Mesh(Mesh &,Geometry * pGh=0,Mesh* pBTh=0,long maxnbv_in=0 ); //copy operator
    … …  
    110110                        BamgVertex* NearestVertex(Icoor1 i,Icoor1 j) ;
    111111                        Triangle* TriangleFindFromCoord(const I2 & ,Icoor2 [3],Triangle *tstart=0) const;
    112                         void ReadMesh(double* index,double* x,double* y,int nods,int nels);
     112                        void ReadMesh(int* index,double* x,double* y,int nods,int nels);
    113113                        void ReadMesh(BamgMesh* bamgmesh, BamgOpts* bamgopts);
    114114                        void WriteMesh(BamgMesh* bamgmesh,BamgOpts* bamgopts);
Note: See TracChangeset for help on using the changeset viewer.