Changeset 13220


Ignore:
Timestamp:
09/02/12 01:29:13 (13 years ago)
Author:
Eric.Larour
Message:

CHG: adjustments to compile the Matlab and Python modules when ADOLC is on.
Had to template MeshPartitionx, because it's used in both modules and the core code.
Had to template Contour too, same reason.
Had to default to ISSMPDouble type for modules such as InterpFromMesh2d* InterpFromGridToMeshx*, etc ..., which
only show up in the Matlab modules.

Location:
issm/trunk-jpl/src
Files:
4 deleted
36 edited

Legend:

Unmodified
Added
Removed

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

    r13216 r13220  
    4646                                        ./classes/objects/Load.h\
    4747                                        ./classes/objects/Contour.h\
    48                                         ./classes/objects/Contour.cpp\
    4948                                        ./classes/objects/OptArgs.h\
    5049                                        ./classes/objects/OptPars.h\
    … …  
    6261                                        ./classes/objects/Node.cpp\
    6362                                        ./classes/objects/Segment.h\
    64                                         ./classes/objects/Segment.cpp\
    6563                                        ./classes/objects/Vertex.h\
    6664                                        ./classes/objects/Vertex.cpp\
    … …  
    311309                                        ./modules/UpdateConstraintsx/UpdateConstraintsx.h\
    312310                                        ./modules/UpdateConstraintsx/UpdateConstraintsx.cpp\
    313                                         ./modules/MeshPartitionx/MeshPartitionx.cpp\
    314311                                        ./modules/MeshPartitionx/MeshPartitionx.h\
    315312                                        ./modules/Reducevectorgtofx/Reducevectorgtofx.cpp\
    … …  
    844841                        ./shared/Exp/IsInPoly.cpp\
    845842                        ./shared/Exp/IsInPolySerial.cpp\
    846                         ./shared/Exp/DomainOutlineRead.cpp\
    847843                        ./shared/Exp/DomainOutlineWrite.cpp\
    848844                        ./shared/TriMesh/trimesh.h\
    … …  
    10141010libISSMModules_a_SOURCES = $(module_sources)
    10151011libISSMModules_a_SOURCES += $(bamg_sources)
     1012if KRIGING
    10161013libISSMModules_a_SOURCES += $(kriging_sources)
     1014endif
     1015if KML
    10171016libISSMModules_a_SOURCES += $(kml_sources)
     1017endif
    10181018libISSMModules_a_CXXFLAGS = $(ALLCXXFLAGS)
    10191019libISSMModules_a_LIBADD = ./libISSMCore.a

  • issm/trunk-jpl/src/c/classes/matrix/Matrix.h

    r13216 r13220  
    7878                /*FUNCTION Matrix(int M,int N,IssmDouble sparsity,int in_type){{{*/
    7979                #ifdef _HAVE_PETSC_
    80                 Matrix(int M,int N,IssmDouble sparsity,int in_type=PetscMatType){
     80                Matrix(int M,int N,double sparsity,int in_type=PetscMatType){
    8181                #else
    82                 Matrix(int M,int N,IssmDouble sparsity,int in_type=SeqMatType){
     82                Matrix(int M,int N,double sparsity,int in_type=SeqMatType){
    8383                #endif
    8484
    … …  
    105105                /*FUNCTION Matrix(IssmDouble* serial_mat, int M,int N,IssmDouble sparsity,int in_type){{{*/
    106106                #ifdef _HAVE_PETSC_
    107                 Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity,int in_type=PetscMatType){
     107                Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity,int in_type=PetscMatType){
    108108                #else
    109                 Matrix(IssmDouble* serial_mat,int M,int N,IssmDouble sparsity,int in_type=SeqMatType){
     109                Matrix(IssmPDouble* serial_mat,int M,int N,IssmPDouble sparsity,int in_type=SeqMatType){
    110110                #endif
    111111

  • issm/trunk-jpl/src/c/classes/objects/Contour.h

    r12832 r13220  
    88/*Headers:*/
    99/*{{{*/
     10#include "../../include/include.h"
     11#include "../../shared/Exceptions/exceptions.h"
     12#include "../../shared/Alloc/xNewDelete.h"
     13#include "../../shared/MemOps/xMemCpy.h"
     14#include "../../io/io.h"
     15#include "../../EnumDefinitions/EnumDefinitions.h"
    1016#include "./Object.h"
    11 #include "../../shared/Exceptions/exceptions.h"
    12 #include "../../toolkits/toolkits.h"
    13 #include "../../include/include.h"
    1417/*}}}*/
    1518
     19template <class doubletype>
    1620class Contour: public Object{
    1721
    … …  
    2024                int     id;
    2125                int       nods;  //number of vertices in the contour
    22                 IssmDouble* x;
    23                 IssmDouble* y;
     26                doubletype* x;
     27                doubletype* y;
    2428                bool    closed; //is this contour closed?
    2529
    26                 /*Contour constructors, destructors {{{*/
    27                 Contour();
    28                 Contour(int id, int nods, IssmDouble* x, IssmDouble* y,bool closed);
    29                 ~Contour();
     30                /*Contour constructors, destructors :*/
     31                /*FUNCTION Contour() default constructor {{{*/
     32                Contour(){
     33                        this->id=0;
     34                        this->nods=0;
     35                        this->x=NULL;
     36                        this->y=NULL;
     37                        this->closed=false;
     38                }
    3039                /*}}}*/
    31                 /*Object virtual functions{{{*/
    32                 void  Echo(void);
    33                 void  DeepEcho(void);
    34                 int   Id(void);
    35                 int   MyRank(void);
    36                 int   ObjectEnum(void);
    37                 Object* copy(void);
     40                /*FUNCTION Contour(int pid, int nods, doubletype* x, doubletype* y,bool closed) {{{*/
     41                Contour(int pid,int pnods, doubletype* px, doubletype* py,bool pclosed){
     42
     43                        this->id=pid;
     44                        this->nods=pnods;
     45                        this->closed=pclosed;
     46                        if(nods){
     47                                this->x=xNew<doubletype>(nods);
     48                                xMemCpy<doubletype>(this->x,px,nods);
     49                                this->y=xNew<doubletype>(nods);
     50                                xMemCpy<doubletype>(this->y,py,nods);
     51                        }
     52                }
     53                /*}}}*/
     54                /*FUNCTION Contour() default constructor {{{*/
     55                ~Contour(){
     56                        xDelete<doubletype>(this->x);
     57                        xDelete<doubletype>(this->y);
     58                }
    3859                /*}}}*/
    3960
    4061
     62                /*Object virtual function resolutoin: */
     63                /*FUNCTION Echo(){{{*/
     64                void Echo(void){
     65
     66                        int i;
     67
     68                        _printLine_("Contour: " << id);
     69                        _printLine_("   nods: " << nods);
     70                        _printLine_("   closed: " << (closed?"true":"false"));
     71                        if(nods){
     72                                _printLine_("   x,y:");
     73                                for(i=0;i<nods;i++){
     74                                        _printLine_(i << ": " << x[i] << "|" << y[i]);
     75                                }
     76                        }
     77                }
     78                /*}}}*/
     79                /*FUNCTION DeepEcho(){{{*/
     80                void DeepEcho(void){
     81                        this->Echo();
     82                }
     83                /*}}}*/
     84                /*FUNCTION Id(){{{*/
     85                int Id(void){
     86                        return id;
     87                }
     88                /*}}}*/
     89                /*FUNCTION MyRank{{{*/
     90                int    MyRank(void){
     91                        extern int my_rank;
     92
     93                        return my_rank;
     94                }
     95                /*}}}*/
     96                /*FUNCTION ObjectEnum{{{*/
     97                int ObjectEnum(void){
     98
     99                        return ContourEnum;
     100
     101                }
     102                /*}}}*/
     103                /*FUNCTION copy {{{*/
     104                Object* copy() {
     105
     106                        return new Contour(*this);
     107
     108                }
     109                /*}}}*/
    41110};
    42111

  • issm/trunk-jpl/src/c/classes/objects/Segment.h

    r12821 r13220  
    1111/*}}}*/
    1212
     13template <class doubletype>
    1314class Segment: public Object{
    1415
    1516        public:
    1617                int eid;
    17                 IssmDouble x1;
    18                 IssmDouble y1;
    19                 IssmDouble x2;
    20                 IssmDouble y2;
     18                doubletype x1;
     19                doubletype y1;
     20                doubletype x2;
     21                doubletype y2;
    2122
    22                 /*Segment constructors, destructors {{{*/
    23                 Segment();
    24                 Segment(int eid,IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2);
    25                 ~Segment();
     23
     24                /*Segment constructors, destructors :*/
     25                /*FUNCTION Segment() default constructor {{{*/
     26                Segment(){
     27                        this->eid=UNDEF;
     28                        this->x1=UNDEF;
     29                        this->y1=UNDEF;
     30                        this->x2=UNDEF;
     31                        this->y2=UNDEF;
     32                }
    2633                /*}}}*/
    27                 /*Object virtual functions definitions:{{{ */
    28                 void  Echo();
    29                 void  DeepEcho();
    30                 int   Id();
    31                 int   MyRank();
    32                 int   ObjectEnum();
    33                 Object* copy();
     34                /*FUNCTION Segment(int eid, doubletype x1,doubletype y1,doubletype x2, doubletype y2){{{*/
     35                Segment(int segment_eid, doubletype segment_x1,doubletype segment_y1,doubletype segment_x2, doubletype segment_y2){
     36
     37                        this->eid=segment_eid;
     38                        this->x1=segment_x1;
     39                        this->y1=segment_y1;
     40                        this->x2=segment_x2;
     41                        this->y2=segment_y2;
     42
     43                }
    3444                /*}}}*/
     45                /*FUNCTION ~Segment(){{{*/
     46                ~Segment(){
     47                }
     48                /*}}}*/
     49
     50                /*Object virtual functions definitions:*/
     51                /*FUNCTION Echo{{{*/
     52                void Echo(void){
     53
     54                        _printLine_("Segment:");
     55                        _printLine_("   eid: " << eid);
     56                        _printLine_("   node 1: " << this->x1 << "|" << this->y1);
     57                        _printLine_("   node 2: " << this->x2 << "|" << this->y2);
     58
     59                }
     60                /*}}}*/
     61                /*FUNCTION DeepEcho{{{*/
     62                void DeepEcho(void){
     63                        this->Echo();
     64                }
     65                /*}}}*/
     66                /*FUNCTION Id{{{*/
     67                int    Id(void){ return eid; }
     68                /*}}}*/
     69                /*FUNCTION MyRank{{{*/
     70                int    MyRank(void){
     71                        extern int my_rank;
     72
     73                        return my_rank;
     74                }
     75                /*}}}*/
     76                /*FUNCTION ObjectEnum{{{*/
     77                int ObjectEnum(void){
     78
     79                        return SegmentEnum;
     80
     81                }
     82                /*}}}*/
     83                /*FUNCTION copy {{{*/
     84                Object* copy() {
     85                        return new Segment(this->eid,this->x1,this->y1,this->x2,this->y2);
     86
     87                }
     88
     89                /*}}}*/
     90
    3591};
    3692

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

    r12832 r13220  
    99#include "../../classes/objects/objects.h"
    1010               
    11 void AddExternalResultx( DataSet* results, int enumtype, double value){
     11void AddExternalResultx( DataSet* results, int enumtype, IssmDouble value){
    1212        /* Add new result in into results*/
    1313        results->AddObject(new DoubleExternalResult(results->Size()+1,enumtype,value,1,0));
    1414}
    15 void AddExternalResultx( DataSet* results, int enumtype, double* value, int nraws){
     15void AddExternalResultx( DataSet* results, int enumtype, IssmDouble* value, int nraws){
    1616        /* Add new result in into results*/
    1717        results->AddObject(new DoubleVecExternalResult(results->Size()+1,enumtype,value,nraws,1,0));

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

    r4554 r13220  
    99
    1010/* local prototypes: */
    11 void AddExternalResultx(DataSet* results, int type, double value);
    12 void AddExternalResultx(DataSet* results, int type, double* value, int nraws);
     11void AddExternalResultx(DataSet* results, int type, IssmDouble value);
     12void AddExternalResultx(DataSet* results, int type, IssmDouble* value, int nraws);
    1313
    1414#endif  /* _ADDEXTERNALRESULT_H */

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

    r13216 r13220  
    1111#include "./ContourToMeshx.h"
    1212
    13 int ContourToMeshx( Vector<IssmDouble>** pin_nod,Vector<IssmDouble>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue) {
     13int ContourToMeshx( Vector<double>** pin_nod,Vector<double>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue) {
    1414
    1515        int noerr=1;
    … …  
    3030
    3131        /*output: */
    32         Vector<IssmDouble>* in_nod=NULL;
    33         Vector<IssmDouble>* in_elem=NULL;
     32        Vector<double>* in_nod=NULL;
     33        Vector<double>* in_elem=NULL;
    3434
    35         in_nod=new Vector<IssmDouble>(nods);
    36         in_elem=new Vector<IssmDouble>(nel);
     35        in_nod=new Vector<double>(nods);
     36        in_elem=new Vector<double>(nel);
    3737
    3838        /*initialize thread parameters: */

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

    r13216 r13220  
    1616        int nods;
    1717        int edgevalue;
    18         Vector<IssmDouble>* in_nod;
     18        Vector<double>* in_nod;
    1919        double* x;
    2020        double* y;
    … …  
    2424
    2525/* local prototypes: */
    26 int ContourToMeshx( Vector<IssmDouble>** pin_nods,Vector<IssmDouble>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue);
     26int ContourToMeshx( Vector<double>** pin_nods,Vector<double>** pin_elem, double* index, double* x, double* y,DataSet* contours,char* interptype,int nel,int nods, int edgevalue);
    2727
    2828void* ContourToMeshxt(void* vContourToMeshxThreadStruct);

  • issm/trunk-jpl/src/c/modules/ContourToMeshx/ContourToMeshxt.cpp

    r13216 r13220  
    3333        double* x=NULL;
    3434        double* y=NULL;
    35         Vector<IssmDouble>* in_nod=NULL;
     35        Vector<double>* in_nod=NULL;
    3636
    3737
    … …  
    5555        /*Loop through all contours: */
    5656        for (i=0;i<contours->Size();i++){
    57                 Contour* contour=(Contour*)contours->GetObjectByOffset(i);
     57                Contour<double>* contour=(Contour<double>*)contours->GetObjectByOffset(i);
    5858                IsInPoly(in_nod,contour->x,contour->y,contour->nods,x,y,i0,i1,edgevalue);
    5959        }

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

    r13216 r13220  
    44#include "./ContourToNodesx.h"
    55
    6 int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, Contour** contours,int numcontours,int edgevalue){
     6int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, Contour<IssmPDouble>** contours,int numcontours,int edgevalue){
    77
    88        int i;
    … …  
    1010
    1111        /*Contour:*/
    12         Contour* contouri=NULL;
     12        Contour<IssmPDouble>* contouri=NULL;
    1313        int      numnodes;
    1414        double*  xc=NULL;
    … …  
    1717
    1818        /*output: */
    19         Vector<IssmDouble>* flags=NULL;
     19        Vector<IssmPDouble>* flags=NULL;
    2020
    21         flags=new Vector<IssmDouble>(nods);
     21        flags=new Vector<IssmPDouble>(nods);
    2222
    2323        /*Loop through all contours: */
    … …  
    3939}
    4040
    41 int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue){
     41int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue){
    4242
    4343        int i;
    … …  
    4545
    4646        /*Contour:*/
    47         Contour* contouri=NULL;
     47        Contour<IssmPDouble>* contouri=NULL;
    4848        int      numnodes;
    4949        double*  xc=NULL;
    … …  
    5252
    5353        /*output: */
    54         Vector<IssmDouble>* flags=NULL;
     54        Vector<IssmPDouble>* flags=NULL;
    5555
    56         flags=new Vector<IssmDouble>(nods);
     56        flags=new Vector<IssmPDouble>(nods);
    5757
    5858        /*Loop through all contours: */
    5959        if(contours){
    6060                for (i=0;i<contours->Size();i++){
    61                         Contour* contour=(Contour*)contours->GetObjectByOffset(i);
     61                        Contour<IssmPDouble>* contour=(Contour<IssmPDouble>*)contours->GetObjectByOffset(i);
    6262                        IsInPoly(flags,contour->x,contour->y,contour->nods,x,y,0,nods,edgevalue);
    6363                }

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

    r13216 r13220  
    1111
    1212/* local prototypes: */
    13 int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, Contour** contours,int numcontours,int edgevalue);
    14 int ContourToNodesx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue);
     13int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, Contour<IssmPDouble>** contours,int numcontours,int edgevalue);
     14int ContourToNodesx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, DataSet* contours, int edgevalue);
    1515
    1616#endif /* _CONTOURTONODESX_H */

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

    r13216 r13220  
    1717
    1818/*InterpFromGridToMeshx{{{*/
    19 int InterpFromGridToMeshx( Vector<IssmDouble>** pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value, int interpenum){
     19int InterpFromGridToMeshx( Vector<IssmPDouble>** pdata_mesh,double* x_in, int x_rows, double* y_in, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods,double default_value, int interpenum){
    2020
    2121        /*output: */
    22         Vector<IssmDouble>* data_mesh=NULL;
     22        Vector<IssmPDouble>* data_mesh=NULL;
    2323       
    2424        /*Intermediary*/
    … …  
    4747
    4848        /*Allocate output vector: */
    49         data_mesh=new Vector<IssmDouble>(nods);
     49        data_mesh=new Vector<IssmPDouble>(nods);
    5050
    5151        /*Find out what kind of coordinates (x_in,y_in) have been given is input*/
    … …  
    127127        double *y             = gate->y;
    128128        int     nods          = gate->nods;
    129         Vector<IssmDouble>*data_mesh     = gate->data_mesh;
     129        Vector<IssmPDouble>*data_mesh     = gate->data_mesh;
    130130        double *data          = gate->data;
    131131        double  default_value = gate->default_value;
    … …  
    181181                                        return NULL; /*WARNING: no error because it would blow up the multithreading!*/
    182182                        }
    183                         if(xIsNan<IssmDouble>(data_value)) data_value=default_value;
     183                        if(xIsNan<IssmPDouble>(data_value)) data_value=default_value;
    184184                }
    185185                else{

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

    r13216 r13220  
    2424        double* x_mesh;
    2525        double* y_mesh;
    26         Vector<IssmDouble>*     data_mesh;
     26        Vector<IssmPDouble>*     data_mesh;
    2727} InterpFromGridToMeshxThreadStruct;
    2828
    29 int    InterpFromGridToMeshx( Vector<IssmDouble>** pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value, int interpenum=BilinearInterpEnum);
     29int    InterpFromGridToMeshx( Vector<IssmPDouble>** pdata_mesh,double* x, int x_rows, double* y, int y_rows, double* data, int M, int N, double* x_mesh, double* y_mesh, int nods, double default_value, int interpenum=BilinearInterpEnum);
    3030void*  InterpFromGridToMeshxt(void* vInterpFromGridToMeshxThreadStruct);
    3131bool   findindices(int* pn,int* pm,double* x,int x_rows, double* y,int y_rows, double xgrid,double ygrid);

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

    r13216 r13220  
    1010#include "../modules.h"
    1111
    12 int InterpFromMesh2dx( Vector<IssmDouble>** pdata_prime,double* 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,
    13                 double* default_values,int num_default_values,Contour** contours,int numcontours){
     12int InterpFromMesh2dx( Vector<IssmPDouble>** pdata_prime,double* 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,
     13                double* default_values,int num_default_values,Contour<IssmPDouble>** contours,int numcontours){
    1414       
    1515        /*Output*/
    16         Vector<IssmDouble>* data_prime=NULL;
     16        Vector<IssmPDouble>* data_prime=NULL;
    1717
    1818        /*Intermediary*/
    … …  
    2727
    2828        /*contours: */
    29         Vector<IssmDouble>*    vec_incontour=NULL;
     29        Vector<IssmPDouble>*    vec_incontour=NULL;
    3030        double*    incontour=NULL;
    3131
    … …  
    7171
    7272        /*Initialize output*/
    73         data_prime=new Vector<IssmDouble>(nods_prime,false,SeqVecType);
     73        data_prime=new Vector<IssmPDouble>(nods_prime,false,SeqVecType);
    7474        if(num_default_values){
    7575                if(num_default_values==1)for (i=0;i<nods_prime;i++) data_prime->SetValue(i,default_values[0],INS_VAL);

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

    r13216 r13220  
    2323        double ymin,ymax;
    2424        int    nods_prime;
    25         Vector<IssmDouble>*    data_prime;
     25        Vector<IssmPDouble>*    data_prime;
    2626        double* x_prime;
    2727        double* y_prime;
    … …  
    3333} InterpFromMesh2dxThreadStruct;
    3434
    35 int InterpFromMesh2dx( Vector<IssmDouble>** pdata_prime,double* 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,
    36                 double* default_values,int num_default_values,Contour** contours,int numcontours);
     35int InterpFromMesh2dx( Vector<IssmPDouble>** pdata_prime,double* 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,
     36                double* default_values,int num_default_values,Contour<IssmPDouble>** contours,int numcontours);
    3737
    3838void* InterpFromMesh2dxt(void* vInterpFromMesh2dxThreadStruct);

  • issm/trunk-jpl/src/c/modules/InterpFromMesh2dx/InterpFromMesh2dxt.cpp

    r13216 r13220  
    3232        double  ymax               = gate->ymax;
    3333        int     nods_prime         = gate->nods_prime;
    34         Vector<IssmDouble>* data_prime         = gate->data_prime;
     34        Vector<IssmPDouble>* data_prime         = gate->data_prime;
    3535        double *x_prime            = gate->x_prime;
    3636        double *y_prime            = gate->y_prime;
    … …  
    8787                                                data_value=data[i];
    8888                                        }
    89                                         if (xIsNan<IssmDouble>(data_value)){
     89                                        if (xIsNan<IssmPDouble>(data_value)){
    9090                                                if(num_default_values==1) data_value=default_values[0];
    9191                                                else data_value=default_values[j];

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

    r13216 r13220  
    77#include "../../include/include.h"
    88
    9 int InterpFromMeshToMesh3dx( Vector<IssmDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value) {
     9int InterpFromMeshToMesh3dx( Vector<IssmPDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value) {
    1010
    1111        /*Output*/
    12         Vector<IssmDouble>* data_prime=NULL;
     12        Vector<IssmPDouble>* data_prime=NULL;
    1313
    1414        /*Intermediary*/
    … …  
    5454
    5555        /*Initialize output*/
    56         data_prime=new Vector<IssmDouble>(nods_prime);
     56        data_prime=new Vector<IssmPDouble>(nods_prime);
    5757        for (i=0;i<nods_prime;i++) data_prime->SetValue(i,default_value,INS_VAL);
    5858
    … …  
    124124                                                data_value=data[i];
    125125                                        }
    126                                         if (xIsNan<IssmDouble>(data_value)) data_value=default_value;
     126                                        if (xIsNan<IssmPDouble>(data_value)) data_value=default_value;
    127127
    128128                                        /*insert value and go to the next point*/

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

    r13216 r13220  
    99#include "../../classes/objects/objects.h"
    1010
    11 int InterpFromMeshToMesh3dx( Vector<IssmDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value);
     11int InterpFromMeshToMesh3dx( Vector<IssmPDouble>** pdata_prime,double* index_data, double* x_data, double* y_data, double* z_data, int nods_data,int nels_data, double* data, int data_length, double* x_prime, double* y_prime, double* z_prime, int nods_prime,double default_value);
    1212
    1313#endif /* _INTERPFROMMESHTOMESH3DX_H */

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

    r12470 r13220  
    77
    88#include "../../include/include.h"
     9#include "../../EnumDefinitions/EnumDefinitions.h"
     10
    911
    1012/* local prototypes: */
    11 int MeshPartitionx(int** pepart, int** pnpart, int numberofelements,int numberofnodes,IssmDouble* elements,
    12                 int numberofelements2d,int numberofnodes2d,IssmDouble* elements2d,int numlayers,int elements_width, int dim,int numareas);
     13template <class doubletype>
     14int MeshPartitionx(int** pepart, int** pnpart, int numberofelements,int numberofnodes,doubletype* elements,
     15                int numberofelements2d,int numberofnodes2d,doubletype* elements2d,int numlayers,int elements_width, int dim,int num_procs){
     16
     17        int noerr=1;
     18        int i,j;
     19
     20        /*Metis partitioning: */
     21        int* epart=NULL;
     22        int* npart=NULL;
     23        int* index=NULL;
     24
     25        int* epart2d=NULL;
     26        int* npart2d=NULL;
     27        int* index2d=NULL;
     28        int  count=0;
     29
     30        int  etype=1; //tria mesh see metis/Programs/Io.c
     31        int  etype2d=1; //tria mesh see metis/Programs/Io.c
     32        int  numflag=0;
     33        int  edgecut=1;
     34
     35        if(dim==2){
     36                epart=xNew<int>(numberofelements);
     37                npart=xNew<int>(numberofnodes);
     38                index=xNew<int>(elements_width*numberofelements);
     39                for (i=0;i<numberofelements;i++){
     40                        for (j=0;j<elements_width;j++){
     41                                *(index+elements_width*i+j)=reCast<int>(*(elements+elements_width*i+j))-1; //-1 for C indexing in Metis
     42                        }
     43                }
     44
     45                /*Partition using Metis:*/
     46                if (num_procs>1){
     47                        #ifdef _HAVE_METIS_
     48                        METIS_PartMeshNodalPatch(&numberofelements,&numberofnodes, index, &etype, &numflag, &num_procs, &edgecut, epart, npart);
     49                        #endif
     50                }
     51                else if (num_procs==1){
     52                        /*METIS does not know how to deal with one cpu only!*/
     53                        for (i=0;i<numberofelements;i++) epart[i]=0;
     54                        for (i=0;i<numberofnodes;i++)    npart[i]=0;
     55                }
     56                else _error_("At least one processor is required");
     57        }
     58        else{
     59                /*We have a 3d mesh, made of a regularly extruded 2d mesh. We first partition the 2d mesh, then we extrude the partition: */
     60
     61                /*First build concatenated 2d mesh  from 2d_coll and 2d_noncoll: */
     62                epart2d=xNew<int>(numberofelements2d);
     63                npart2d=xNew<int>(numberofnodes2d);
     64                index2d=xNew<int>(3*numberofelements2d);
     65
     66                for (i=0;i<numberofelements2d;i++){
     67                        for (j=0;j<3;j++){
     68                                *(index2d+3*i+j)=reCast<int>(*(elements2d+3*i+j))-1; //-1 for C indexing in Metis
     69                        }
     70                }
     71
     72                /*Partition using Metis:*/
     73                if (num_procs>1){
     74                        #ifdef _HAVE_METIS_
     75                        METIS_PartMeshNodalPatch(&numberofelements2d,&numberofnodes2d, index2d, &etype2d, &numflag, &num_procs, &edgecut, epart2d, npart2d);
     76                        #endif
     77                }
     78                else if (num_procs==1){
     79                        /*METIS does not know how to deal with one cpu only!*/
     80                        for (i=0;i<numberofelements2d;i++) epart2d[i]=0;
     81                        for (i=0;i<numberofnodes2d;i++)    npart2d[i]=0;
     82                }
     83                else _error_("At least one processor is required");
     84
     85                /*Extrude epart2d to epart, using numlayers: */
     86                epart=xNew<int>(numberofelements);
     87               
     88                count=0;
     89                for(i=0;i<(numlayers-1);i++){
     90                        for(j=0;j<numberofelements2d;j++){
     91                                epart[count]=epart2d[j];
     92                                count++;
     93                        }
     94                }
     95
     96                /*Extrude npart2d to npart, using numlayers: */
     97                npart=xNew<int>(numberofnodes);
     98               
     99                count=0;
     100                for(i=0;i<(numlayers);i++){
     101                        for(j=0;j<numberofnodes2d;j++){
     102                                npart[count]=npart2d[j];
     103                                count++;
     104                        }
     105                }
     106        }
    13107       
     108        /*Assign output pointer:*/
     109        *pepart=epart;
     110        *pnpart=npart;
     111
     112        /*Free ressources: */
     113        xDelete<int>(index);
     114        xDelete<int>(epart2d);
     115        xDelete<int>(npart2d);
     116        xDelete<int>(index2d);
     117        return noerr;
     118}       
    14119#endif /* _MESHPARTITIONX_H */
    15120

  • issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/ElementSegment.cpp

    r13056 r13220  
    5555                        yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
    5656
    57                         segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
     57                        segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
    5858                }
    5959                else{
    … …  
    8383                yfinal[1]=ysegment[0]+coord2*(ysegment[1]-ysegment[0]);
    8484
    85                 segments_dataset->AddObject(new  Segment(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
     85                segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
    8686        }
    8787        else{
    8888                /*No interesections, but the segment might be entirely inside this triangle!: */
    8989                if ( (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) && (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])) ){
    90                         segments_dataset->AddObject(new  Segment(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
     90                        segments_dataset->AddObject(new  Segment<double>(el+1,xsegment[0],ysegment[0],xsegment[1],ysegment[1]));
    9191                }
    9292        }

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

    r12450 r13220  
    44#include "./MeshProfileIntersectionx.h"
    55
    6 void MeshProfileIntersectionx( double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods,  Contour** contours,int numcontours){
     6void MeshProfileIntersectionx( double** psegments, int* pnumsegs, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours){
    77
    88        int i,j,k;
    … …  
    1010
    1111        /*Contour:*/
    12         Contour* contouri=NULL;
     12        Contour<IssmPDouble>* contouri=NULL;
    1313        int      numnodes;
    1414        double*  xc=NULL;

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

    r12832 r13220  
    1111
    1212/* local prototypes: */
    13 void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour** contours,int numcontours);
     13void MeshProfileIntersectionx( double** psegments, int* pnumseg, int* index, double* x, double* y, int nel, int nods,  Contour<IssmPDouble>** contours,int numcontours);
    1414void MeshSegmentsIntersection(double** psegments, int* pnumsegs,int* index, double* x, double* y, int nel, int nods, double* xc, double* yc, int numnodes);
    1515void ElementSegmentsIntersection(DataSet* segments_dataset,int el, double* xnodes,double* ynodes,double* xc,double* yc,int numnodes);

  • issm/trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshSegmentsIntersection.cpp

    r12450 r13220  
    1010        /*output: */
    1111        double*  segments=NULL;
    12         Segment* segment=NULL;
     12        Segment<double>* segment=NULL;
    1313        int     numsegs;
    1414       
    … …  
    3434        segments=xNew<double>(5*numsegs);
    3535        for(i=0;i<numsegs;i++){
    36                 Segment* segment=(Segment*)segments_dataset->GetObjectByOffset(i);
     36                Segment<double>* segment=(Segment<double>*)segments_dataset->GetObjectByOffset(i);
    3737               
    3838                /*x1,y1,x2,y2 then element_id: */

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

    r13216 r13220  
    44#include "./PointCloudFindNeighborsx.h"
    55
    6 int PointCloudFindNeighborsx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread){
     6int PointCloudFindNeighborsx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread){
    77
    88        /*output: */
    9         Vector<IssmDouble>* flags=NULL;
    10         flags=new Vector<IssmDouble>(nods);
     9        Vector<IssmPDouble>* flags=NULL;
     10        flags=new Vector<IssmPDouble>(nods);
    1111
    1212        /*threading: */

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

    r13216 r13220  
    1111
    1212/* local prototypes: */
    13 int PointCloudFindNeighborsx( Vector<IssmDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread);
     13int PointCloudFindNeighborsx( Vector<IssmPDouble>** pflags,double* x, double* y, int nods, double mindistance,double multithread);
    1414
    1515/*threading: */
    … …  
    2020        int nods;
    2121        double mindistance;
    22         Vector<IssmDouble>* flags;
     22        Vector<IssmPDouble>* flags;
    2323
    2424

  • issm/trunk-jpl/src/c/modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp

    r13216 r13220  
    1717        int     nods;
    1818        double  mindistance;
    19         Vector<IssmDouble>*     flags;
     19        Vector<IssmPDouble>*     flags;
    2020
    2121        /*recover handle and gate: */

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

    r13216 r13220  
    1515        int configuration_type;
    1616        int ssize;
    17         double* yg_serial=NULL;
     17        IssmDouble* yg_serial=NULL;
    1818
    1919        /*first figure out ssize: */
    … …  
    5151
    5252        /*Free ressources:*/
    53         xDelete<double>(yg_serial);
     53        xDelete<IssmDouble>(yg_serial);
    5454
    5555        /*Assign output pointers:*/

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

    r13216 r13220  
    2020/*}}}*/
    2121
    22 void TriMeshx(Matrix<IssmDouble>** pindex,Vector<IssmDouble>** px,Vector<IssmDouble>** py,Matrix<IssmDouble>** psegments,Vector<IssmDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area){
     22void TriMeshx(Matrix<IssmPDouble>** pindex,Vector<IssmPDouble>** px,Vector<IssmPDouble>** py,Matrix<IssmPDouble>** psegments,Vector<IssmPDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area){
    2323
    2424        /*indexing: */
    … …  
    2626
    2727        /*output: */
    28         double* index=NULL;
    29         Matrix<IssmDouble>* index_matrix=NULL;
     28        IssmPDouble* index=NULL;
     29        Matrix<IssmPDouble>* index_matrix=NULL;
    3030        double* x=NULL;
    3131        double* y=NULL;
    3232        double* segments=NULL;
    33         Matrix<IssmDouble>* segments_matrix=NULL;
     33        Matrix<IssmPDouble>* segments_matrix=NULL;
    3434        double* segmentmarkerlist=NULL;
    3535
    3636        /*intermediary: */
    3737        int      counter,counter2,backcounter;
    38         Contour* contour=NULL;
     38        Contour<IssmPDouble>* contour=NULL;
    3939
    4040        /* Triangle structures needed to call Triangle library routines: */
    … …  
    4545        in.numberofpoints=0;
    4646        for (i=0;i<domain->Size();i++){
    47                 contour=(Contour*)domain->GetObjectByOffset(i);
     47                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    4848                in.numberofpoints+=contour->nods-1;
    4949        }
    5050        for (i=0;i<rifts->Size();i++){
    51                 contour=(Contour*)rifts->GetObjectByOffset(i);
     51                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    5252                in.numberofpoints+=contour->nods;
    5353        }
    … …  
    6161        counter=0;
    6262        for (i=0;i<domain->Size();i++){
    63                 contour=(Contour*)domain->GetObjectByOffset(i);
     63                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    6464                for (j=0;j<contour->nods-1;j++){
    6565                        in.pointlist[2*counter+0]=contour->x[j];
    … …  
    6969        }
    7070        for (i=0;i<rifts->Size();i++){
    71                 contour=(Contour*)rifts->GetObjectByOffset(i);
     71                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    7272                for (j=0;j<contour->nods;j++){
    7373                        in.pointlist[2*counter+0]=contour->x[j];
    … …  
    8888        in.numberofsegments=0;
    8989        for (i=0;i<domain->Size();i++){
    90                 contour=(Contour*)domain->GetObjectByOffset(i);
     90                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    9191                in.numberofsegments+=contour->nods-1;
    9292        }
    9393        for(i=0;i<rifts->Size();i++){
    94                 contour=(Contour*)rifts->GetObjectByOffset(i);
     94                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    9595                /*for rifts, we have one less segment as we have vertices*/
    9696                in.numberofsegments+=contour->nods-1;
    … …  
    102102        backcounter=0;
    103103        for (i=0;i<domain->Size();i++){
    104                 contour=(Contour*)domain->GetObjectByOffset(i);
     104                contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i);
    105105                for (j=0;j<contour->nods-2;j++){
    106106                        in.segmentlist[2*counter+0]=counter;
    … …  
    118118        counter2=counter;
    119119        for (i=0;i<rifts->Size();i++){
    120                 contour=(Contour*)rifts->GetObjectByOffset(i);
     120                contour=(Contour<IssmPDouble>*)rifts->GetObjectByOffset(i);
    121121                for (j=0;j<(contour->nods-1);j++){
    122122                        in.segmentlist[2*counter2+0]=counter;
    … …  
    137137                in.holelist = xNew<REAL>(in.numberofholes*2);
    138138                for (i=0;i<domain->Size()-1;i++){
    139                         contour=(Contour*)domain->GetObjectByOffset(i+1);
     139                        contour=(Contour<IssmPDouble>*)domain->GetObjectByOffset(i+1);
    140140                        GridInsideHole(&in.holelist[2*i+0],&in.holelist[2*i+1],contour->nods-1,contour->x,contour->y);
    141141                }
    … …  
    192192
    193193        /*Output : */
    194         index_matrix=new Matrix<IssmDouble>(index,out.numberoftriangles,3,1,SeqMatType);
     194        index_matrix=new Matrix<IssmPDouble>(index,out.numberoftriangles,3,1.0,SeqMatType);
    195195        *pindex=index_matrix;
    196196       
    197         segments_matrix=new Matrix<IssmDouble>(segments,out.numberofsegments,3,1,SeqMatType);
     197        segments_matrix=new Matrix<IssmPDouble>(segments,out.numberofsegments,3,1.0,SeqMatType);
    198198        *psegments=segments_matrix;
    199199
    200         *px=new Vector<IssmDouble>(x,out.numberofpoints,SeqMatType);
    201         *py=new Vector<IssmDouble>(y,out.numberofpoints,SeqMatType);
    202         *psegmentmarkerlist=new Vector<IssmDouble>(segmentmarkerlist,out.numberofsegments,SeqMatType);
     200        *px=new Vector<IssmPDouble>(x,out.numberofpoints,SeqMatType);
     201        *py=new Vector<IssmPDouble>(y,out.numberofpoints,SeqMatType);
     202        *psegmentmarkerlist=new Vector<IssmPDouble>(segmentmarkerlist,out.numberofsegments,SeqMatType);
    203203}

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

    r13216 r13220  
    1111
    1212/* local prototypes: */
    13 void TriMeshx(Matrix<IssmDouble>** pindex,Vector<IssmDouble>** px,Vector<IssmDouble>** py,Matrix<IssmDouble>** psegments,Vector<IssmDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area);
     13void TriMeshx(Matrix<IssmPDouble>** pindex,Vector<IssmPDouble>** px,Vector<IssmPDouble>** py,Matrix<IssmPDouble>** psegments,Vector<IssmPDouble>** psegmentmarkerlist,DataSet* domain,DataSet* rifts,double area);
    1414
    1515#endif  /* _TRIMESHX_H */

  • issm/trunk-jpl/src/c/shared/Exp/IsInPoly.cpp

    r13216 r13220  
    1414#endif
    1515
    16 /*IsInPoly {{{*/
    17 int IsInPoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
    18 
    19         int i;
    20         double x0,y0;
    21         double value;
    22         double xmin=xc[0];
    23         double xmax=xc[0];
    24         double ymin=yc[0];
    25         double ymax=yc[0];
    26 
    27         /*Get extrema*/
    28         for (i=1;i<numvertices;i++){
    29                 if(xc[i]<xmin) xmin=xc[i];
    30                 if(xc[i]>xmax) xmax=xc[i];
    31                 if(yc[i]<ymin) ymin=yc[i];
    32                 if(yc[i]>ymax) ymax=yc[i];
    33         }
    34 
    35         /*Go through all vertices of the mesh:*/
    36         for (i=i0;i<i1;i++){
    37 
    38                 //Get current value of value[i] -> do not change it if != 0
    39                 in->GetValue(&value,i);
    40                 if (value){
    41                         /*this vertex already is inside one of the contours, continue*/
    42                         continue;
    43                 }
    44 
    45                 /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
    46                 x0=x[i]; y0=y[i];
    47                 if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
    48                         value=0;
    49                 }
    50                 else{
    51                         value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
    52                 }
    53                 in->SetValue(i,value,INS_VAL);
    54         }
    55          return 1;
    56 }/*}}}*/
    57 /*IsOutsidePoly {{{*/
    58 int IsOutsidePoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
    59 
    60         int i,j;
    61         double x0,y0;
    62         double value;
    63         double xmin=xc[0];
    64         double xmax=xc[0];
    65         double ymin=yc[0];
    66         double ymax=yc[0];
    67 
    68         /*Get extrema*/
    69         for (i=1;i<numvertices;i++){
    70                 if(xc[i]<xmin) xmin=xc[i];
    71                 if(xc[i]>xmax) xmax=xc[i];
    72                 if(yc[i]<ymin) ymin=yc[i];
    73                 if(yc[i]>ymax) ymax=yc[i];
    74         }
    75 
    76         /*Go through all vertices of the mesh:*/
    77         for (i=i0;i<i1;i++){
    78 
    79                 //Get current value of value[i] -> do not change it if != 0
    80                 in->GetValue(&value,i);
    81                 if (value){
    82                         /*this vertex already is inside one of the contours, continue*/
    83                         continue;
    84                 }
    85 
    86                 /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
    87                 x0=x[i]; y0=y[i];
    88                 if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
    89                         value=1;
    90                 }
    91                 else{
    92                         value=1-pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
    93                 }
    94                 in->SetValue(i,value,INS_VAL);
    95         }
    96         return 1;
    97 }/*}}}*/
    9816/*pnpoly{{{*/
    9917int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue) {

  • issm/trunk-jpl/src/c/shared/Exp/exp.h

    r13216 r13220  
    77#define _EXP_H_
    88
    9 #include "../../classes/objects/objects.h"
    10 #include "../../toolkits/toolkits.h"
    11 
    12 int IsInPoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue);
    13 int IsOutsidePoly(Vector<IssmDouble>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue);
     9#include "../../classes/objects/Contour.h"
     10#include "../../shared/Numerics/recast.h"
     11#include "../../Container/Container.h"
     12
    1413int IsInPolySerial(double* in,double* xc,double* yc,int numvertices,double* x,double* y,int nods, int edgevalue);
    1514int DomainOutlineWrite(int nprof,int* profnvertices,double** pprofx,double** pprofy,bool* closed,char* domainname);
    1615int pnpoly(int npol, double *xp, double *yp, double x, double y, int edgevalue);
    1716
    18 int      DomainOutlineRead(int* pnprof,int** pprofnvertices,double*** ppprofx,double*** ppprofy,bool** pclosed,char* domainname);
    19 DataSet* DomainOutlineRead(char* domainname);
     17
     18/*IsInPoly {{{*/
     19template <class doubletype>
     20int IsInPoly(Vector<doubletype>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
     21
     22        int i;
     23        double x0,y0;
     24        doubletype value;
     25        double xmin=xc[0];
     26        double xmax=xc[0];
     27        double ymin=yc[0];
     28        double ymax=yc[0];
     29
     30        /*Get extrema*/
     31        for (i=1;i<numvertices;i++){
     32                if(xc[i]<xmin) xmin=xc[i];
     33                if(xc[i]>xmax) xmax=xc[i];
     34                if(yc[i]<ymin) ymin=yc[i];
     35                if(yc[i]>ymax) ymax=yc[i];
     36        }
     37
     38        /*Go through all vertices of the mesh:*/
     39        for (i=i0;i<i1;i++){
     40
     41                //Get current value of value[i] -> do not change it if != 0
     42                in->GetValue(&value,i);
     43                if (reCast<bool,doubletype>(value)){
     44                        /*this vertex already is inside one of the contours, continue*/
     45                        continue;
     46                }
     47
     48                /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
     49                x0=x[i]; y0=y[i];
     50                if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
     51                        value=0;
     52                }
     53                else{
     54                        value=pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
     55                }
     56                in->SetValue(i,value,INS_VAL);
     57        }
     58         return 1;
     59}/*}}}*/
     60/*IsOutsidePoly {{{*/
     61template <class doubletype>
     62int IsOutsidePoly(Vector<doubletype>* in,double* xc,double* yc,int numvertices,double* x,double* y,int i0,int i1, int edgevalue){
     63
     64        int i,j;
     65        double x0,y0;
     66        doubletype value;
     67        double xmin=xc[0];
     68        double xmax=xc[0];
     69        double ymin=yc[0];
     70        double ymax=yc[0];
     71
     72        /*Get extrema*/
     73        for (i=1;i<numvertices;i++){
     74                if(xc[i]<xmin) xmin=xc[i];
     75                if(xc[i]>xmax) xmax=xc[i];
     76                if(yc[i]<ymin) ymin=yc[i];
     77                if(yc[i]>ymax) ymax=yc[i];
     78        }
     79
     80        /*Go through all vertices of the mesh:*/
     81        for (i=i0;i<i1;i++){
     82
     83                //Get current value of value[i] -> do not change it if != 0
     84                in->GetValue(&value,i);
     85                if (reCast<bool,doubletype>(value)){
     86                        /*this vertex already is inside one of the contours, continue*/
     87                        continue;
     88                }
     89
     90                /*pick up vertex (x[i],y[i]) and figure out if located inside contour (xc,yc)*/
     91                x0=x[i]; y0=y[i];
     92                if(x0<xmin || x0>xmax || y0<ymin || y0>ymax){
     93                        value=1;
     94                }
     95                else{
     96                        value=1-pnpoly(numvertices,xc,yc,x0,y0,edgevalue);
     97                }
     98                in->SetValue(i,value,INS_VAL);
     99        }
     100        return 1;
     101}/*}}}*/
     102/*DomainOutlineRead{{{*/
     103template <class doubletype>
     104int DomainOutlineRead(int* pnprof,int** pprofnvertices,doubletype*** ppprofx,doubletype*** ppprofy,bool** pclosed,char* domainname){
     105       
     106        /*indexing: */
     107        int i,counter;
     108
     109        /*I/O: */
     110        FILE   *fid = NULL;
     111        char    chardummy[256];
     112        double  ddummy;
     113
     114        /*output: */
     115        int      nprof;                //number of profiles in the domainname file
     116        int     *profnvertices = NULL; //array holding the number of vertices for the nprof profiles
     117        doubletype **pprofx        = NULL; //array of profiles x coordinates
     118        doubletype **pprofy        = NULL; //array of profiles y coordinates
     119        bool    *closed        = NULL; //array holding closed flags for the nprof profiles
     120
     121        /*For each profile: */
     122        int     n;
     123        doubletype *x  = NULL;
     124        doubletype *y  = NULL;
     125        double  xi,yi;
     126        bool    cl;
     127
     128        /*open domain outline file for reading: */
     129        if ((fid=fopen(domainname,"r"))==NULL){
     130                _error_("could not find domain file " << domainname);
     131        }
     132
     133        /*Do a first pass through the domainname file, to figure out how many profiles
     134         *we need to read: */
     135        nprof=1;
     136        for(;;){
     137                fscanf(fid,"%256s %256s\n",chardummy,chardummy);
     138                fscanf(fid,"%256s %256s\n",chardummy,chardummy);
     139                fscanf(fid,"%256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy);
     140                fscanf(fid,"%20u %256s\n",&n,chardummy);
     141                fscanf(fid,"%256s %256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
     142                for (i=0;i<n;i++){
     143                        fscanf(fid,"%20lf %20lf\n",&ddummy,&ddummy);
     144                }
     145                /*Ok, we have faked one profile reading, check whether we are at the end of the file, otherwise, keep fake reading next profile:*/
     146                if (feof(fid)){
     147                        break;
     148                }
     149                nprof++;
     150        }
     151       
     152        /*Allocate and initialize all the profiles: */
     153        profnvertices=xNew<int>(nprof);
     154        pprofx=xNew<doubletype*>(nprof);
     155        pprofy=xNew<doubletype*>(nprof);
     156        for (i=0;i<nprof;i++){
     157                pprofx[i]=NULL;
     158                pprofy[i]=NULL;
     159        }
     160        closed=xNew<bool>(nprof);
     161
     162        /*Reaset file pointer to beginning of file: */
     163        fseek(fid,0,SEEK_SET);
     164
     165        /*Start reading profiles: */
     166        for(counter=0;counter<nprof;counter++){
     167
     168                /*Skip header: */
     169                fscanf(fid,"%256s %256s\n",chardummy,chardummy);
     170                fscanf(fid,"%256s %256s\n",chardummy,chardummy);
     171                fscanf(fid,"%256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy);
     172               
     173                /*Get number of profile vertices: */
     174                fscanf(fid,"%20u %256s\n",&n,chardummy);
     175       
     176                /*Skip next line: */
     177                fscanf(fid,"%256s %256s %256s %256s %256s\n",chardummy,chardummy,chardummy,chardummy,chardummy);
     178
     179                /*Allocate vertices: */
     180                x=xNew<doubletype>(n);
     181                y=xNew<doubletype>(n);
     182               
     183                /*Read vertices: */
     184                for (i=0;i<n;i++){
     185                        fscanf(fid,"%20lf %20lf\n",&xi,&yi);
     186                        x[i]=xi; y[i]=yi;
     187                }
     188
     189                /*Now check that we are dealing with open contours: */
     190                cl=false;
     191                if((x[0]==x[n-1]) && (y[0]==y[n-1])){
     192                        cl=true;
     193                }
     194
     195                /*Assign pointers: */
     196                profnvertices[counter]=n;
     197                pprofx[counter]=x;
     198                pprofy[counter]=y;
     199                closed[counter]=cl;
     200        }
     201
     202        /*close domain outline file: */
     203        fclose(fid);
     204
     205        /*Assign output pointers: */
     206        *pnprof=nprof;
     207        *pprofnvertices=profnvertices;
     208        *ppprofx=pprofx;
     209        *ppprofy=pprofy;
     210        if(pclosed)
     211         *pclosed=closed;
     212        else
     213         xDelete<bool>(closed);
     214} /*}}}*/
     215/*DataSet* DomainOutlineRead(char* domainname){{{*/
     216template <class doubletype>
     217DataSet* DomainOutlineRead(char* domainname){
     218
     219        /*intermediary: */
     220        int       nprof;
     221        int      *profnvertices = NULL;
     222        doubletype  **pprofx        = NULL;
     223        doubletype  **pprofy        = NULL;
     224        Contour<doubletype>* contour       = NULL;
     225
     226        /*output: */
     227        DataSet* domain=NULL;
     228
     229        /*If domainname is an empty string, return empty dataset*/
     230        if (strcmp(domainname,"")==0){
     231                nprof=0;
     232        }
     233        else{
     234                DomainOutlineRead(&nprof,&profnvertices,&pprofx, &pprofy, NULL,domainname);
     235        }
     236
     237        /*now create dataset of contours: */
     238        domain=new DataSet(0);
     239
     240        for(int i=0;i<nprof;i++){
     241                domain->AddObject(new Contour<doubletype>(i,profnvertices[i],pprofx[i],pprofy[i],1));
     242        }
     243        return domain;
     244} /*}}}*/
     245
     246
    20247
    21248

  • issm/trunk-jpl/src/modules/ContourToMesh/ContourToMesh.cpp

    r13216 r13220  
    6363        FetchData(&edgevalue,EDGEVALUE);
    6464        FetchData(&contourname,CONTOURNAME);
    65         contours=DomainOutlineRead(contourname);
     65        contours=DomainOutlineRead<double>(contourname);
    6666
    6767        /*Fetch  interptype: */

  • issm/trunk-jpl/src/modules/ContourToNodes/ContourToNodes.cpp

    r13216 r13220  
    4141        mxArray*  matlabstructure=NULL;
    4242        int numcontours;
    43         Contour** contours=NULL;
    44         Contour*  contouri=NULL;
     43        Contour<double>** contours=NULL;
     44        Contour<double>*  contouri=NULL;
    4545
    4646        /*Boot module: */
    … …  
    6868
    6969        numcontours=mxGetNumberOfElements(matlabstructure);
    70         contours=xNew<Contour*>(numcontours);
     70        contours=xNew<Contour<double>*>(numcontours);
    7171        for(i=0;i<numcontours;i++){
    7272                //allocate
    73                 contouri=xNew<Contour>(1);
     73                contouri=xNew<Contour<double> >(1);
    7474                //retrieve dimension of this contour.
    7575                contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));

  • issm/trunk-jpl/src/modules/InterpFromMesh2d/InterpFromMesh2d.cpp

    r13216 r13220  
    5050        //contours
    5151        mxArray*  matlabstructure=NULL;
    52         Contour** contours=NULL;
     52        Contour<double>** contours=NULL;
    5353        int       numcontours;
    54         Contour*  contouri=NULL;
     54        Contour<double>*  contouri=NULL;
    5555        int       i;
    5656
    … …  
    100100                /*contours: */
    101101                numcontours=mxGetNumberOfElements(matlabstructure);
    102                 contours=xNew<Contour*>(numcontours);
     102                contours=xNew<Contour<double> *>(numcontours);
    103103                for(i=0;i<numcontours;i++){
    104104                        //allocate
    105                         contouri=xNew<Contour>(1);
     105                        contouri=xNew<Contour<double> >(1);
    106106                        //retrieve dimension of this contour.
    107107                        contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));

  • issm/trunk-jpl/src/modules/MeshProfileIntersection/MeshProfileIntersection.cpp

    r13038 r13220  
    3434        //contours
    3535        mxArray*  matlabstructure=NULL;
    36         Contour** contours=NULL;
     36        Contour<double>** contours=NULL;
    3737        int       numcontours;
    38         Contour*  contouri=NULL;
     38        Contour<double>*  contouri=NULL;
    3939
    4040        /* output: */
    … …  
    6767        //contours
    6868        numcontours=mxGetNumberOfElements(matlabstructure);
    69         contours=xNew<Contour*>(numcontours);
     69        contours=xNew<Contour<double>*>(numcontours);
    7070        for(i=0;i<numcontours;i++){
    7171                //allocate
    72                 contouri=xNew<Contour>(1);
     72                contouri=xNew<Contour<double> >(1);
    7373                //retrieve dimension of this contour.
    7474                contouri->nods=(int)mxGetScalar(mxGetField(matlabstructure,i,"nods"));

  • issm/trunk-jpl/src/modules/TriMesh/TriMesh.cpp

    r13216 r13220  
    3535
    3636        /*Read domain outline: */
    37         domain = DomainOutlineRead(domainname);
    38         rifts  = DomainOutlineRead(riftsname);
     37        domain = DomainOutlineRead<double>(domainname);
     38        rifts  = DomainOutlineRead<double>(riftsname);
    3939
    4040        /*call x core: */
Note: See TracChangeset for help on using the changeset viewer.