Changeset 304


Ignore:
Timestamp:
05/07/09 12:13:45 (16 years ago)
Author:
Eric.Larour
Message:

New framework to include hutter and stokes formulations

Location:
issm/trunk
Files:
1 added
21 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/BuildNodeSetsx/BuildNodeSetsx.cpp

    r1 r304  
    4040        gsize=nodes->NumberOfDofs();
    4141
    42         /*Now, allocate parallel vectors for all the sets: */
    43         flag_pv_g=NewVec(gsize);
    44         flag_pv_m=NewVec(gsize);
    45         flag_pv_n=NewVec(gsize);
    46         flag_pv_f=NewVec(gsize);
    47         flag_pv_s=NewVec(gsize);
     42        if(gsize){
    4843
    49         /*Now, go through all nodes and have every node plug into
    50          * pv_m, pv_n, pv_f and pv_s, their own set flags, for each dof: */
    51         nodes->FlagNodeSets(flag_pv_g,flag_pv_m,flag_pv_n,flag_pv_f,flag_pv_s);
     44                /*Now, allocate parallel vectors for all the sets: */
     45                flag_pv_g=NewVec(gsize);
     46                flag_pv_m=NewVec(gsize);
     47                flag_pv_n=NewVec(gsize);
     48                flag_pv_f=NewVec(gsize);
     49                flag_pv_s=NewVec(gsize);
    5250
    53         /*Now, every cpu has 4 flag vectors, of size gsize. Create partition vectors (like a pos=find(flag_pv_g) in matlab*/
    54         PartitionSets(&vec_pv_m,&vec_pv_n,flag_pv_g,flag_pv_m,flag_pv_n,gsize); /*! split g set into m and n sets*/
    55         PartitionSets(&vec_pv_f,&vec_pv_s,flag_pv_n,flag_pv_f,flag_pv_s,gsize); /*! split n set into f and s sets*/
     51                /*Now, go through all nodes and have every node plug into
     52                 * pv_m, pv_n, pv_f and pv_s, their own set flags, for each dof: */
     53                nodes->FlagNodeSets(flag_pv_g,flag_pv_m,flag_pv_n,flag_pv_f,flag_pv_s);
    5654
    57         /*Free ressources:*/
    58         VecFree(&flag_pv_g);
    59         VecFree(&flag_pv_m);
    60         VecFree(&flag_pv_n);
    61         VecFree(&flag_pv_f);
    62         VecFree(&flag_pv_s);
     55                /*Now, every cpu has 4 flag vectors, of size gsize. Create partition vectors (like a pos=find(flag_pv_g) in matlab*/
     56                PartitionSets(&vec_pv_m,&vec_pv_n,flag_pv_g,flag_pv_m,flag_pv_n,gsize); /*! split g set into m and n sets*/
     57                PartitionSets(&vec_pv_f,&vec_pv_s,flag_pv_n,flag_pv_f,flag_pv_s,gsize); /*! split n set into f and s sets*/
    6358
    64         /*Serialize vectors: */
    65         VecGetSize(vec_pv_m,&msize);
    66         VecGetSize(vec_pv_n,&nsize);
    67         VecGetSize(vec_pv_f,&fsize);
    68         VecGetSize(vec_pv_s,&ssize);
    69         VecToMPISerial(&pv_m,vec_pv_m);
    70         VecToMPISerial(&pv_n,vec_pv_n);
    71         VecToMPISerial(&pv_f,vec_pv_f);
    72         VecToMPISerial(&pv_s,vec_pv_s);
     59                /*Free ressources:*/
     60                VecFree(&flag_pv_g);
     61                VecFree(&flag_pv_m);
     62                VecFree(&flag_pv_n);
     63                VecFree(&flag_pv_f);
     64                VecFree(&flag_pv_s);
    7365
    74         /*Free ressources:*/
    75         VecFree(&vec_pv_m);
    76         VecFree(&vec_pv_n);
    77         VecFree(&vec_pv_f);
    78         VecFree(&vec_pv_s);
     66                /*Serialize vectors: */
     67                VecGetSize(vec_pv_m,&msize);
     68                VecGetSize(vec_pv_n,&nsize);
     69                VecGetSize(vec_pv_f,&fsize);
     70                VecGetSize(vec_pv_s,&ssize);
     71                VecToMPISerial(&pv_m,vec_pv_m);
     72                VecToMPISerial(&pv_n,vec_pv_n);
     73                VecToMPISerial(&pv_f,vec_pv_f);
     74                VecToMPISerial(&pv_s,vec_pv_s);
     75
     76                /*Free ressources:*/
     77                VecFree(&vec_pv_m);
     78                VecFree(&vec_pv_n);
     79                VecFree(&vec_pv_f);
     80                VecFree(&vec_pv_s);
    7981
    8082
    81         /*Create NodeSets* object: */
    82         nodesets=new NodeSets(pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
     83                /*Create NodeSets* object: */
     84                nodesets=new NodeSets(pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
     85        }
    8386
    8487        /*Assign output pointers: */

  • issm/trunk/src/c/DataSet/DataSet.cpp

    r246 r304  
    150150        vector<Object*>::iterator iterator;
    151151
    152         iterator = find(objects.begin(), objects.end(),object);
    153 
    154         objects.erase(iterator);
     152        if(object){
     153                iterator = find(objects.begin(), objects.end(),object);
     154                objects.erase(iterator);
     155        }
    155156
    156157}
    … …  
    11391140
    11401141}
     1142
     1143
     1144void DataSet::ComputePressure(Vec p_g){
     1145
     1146        vector<Object*>::iterator object;
     1147        Element* element=NULL;
     1148
     1149        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1150               
     1151                if(EnumIsElement((*object)->Enum())){
     1152
     1153                        element=(Element*)(*object);
     1154                        element->ComputePressure(p_g);
     1155                }
     1156        }
     1157
     1158}
     1159
     1160

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

    r246 r304  
    7676                void  VelocityExtrude(Vec ug,double* ug_serial);
    7777                int   DeleteObject(Object* object);
     78                void  ComputePressure(Vec p_g);
    7879
    7980};

  • issm/trunk/src/c/EnumDefinitions/AnalysisTypeAsEnum.cpp

    r128 r304  
    3939                return DiagnosticStokesAnalysisEnum();
    4040        }
     41        else if (strcmp(analysis_type,"diagnostic_hutter")==0){
     42                return DiagnosticHutterAnalysisEnum();
     43        }
     44        else if (strcmp(analysis_type,"surface_slope_compute")==0){
     45                return SurfaceSlopeComputeAnalysisEnum();
     46        }
     47        else if (strcmp(analysis_type,"bed_slope_compute")==0){
     48                return BedSlopeComputeAnalysisEnum();
     49        }
    4150        else throw ErrorException(__FUNCT__," could not recognized analysis type!");
    4251

  • issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp

    r246 r304  
    2525int MeltingAnalysisEnum(void){          return          29; }
    2626int DiagnosticStokesAnalysisEnum(void){ return          30; }
     27int DiagnosticHutterAnalysisEnum(void){ return          31; }
     28int SurfaceSlopeComputeAnalysisEnum(void){ return       32; }
     29int BedSlopeComputeAnalysisEnum(void){ return           33; }
    2730
    2831/*datasets: */

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

    r246 r304  
    3939int MeltingAnalysisEnum(void);
    4040int DiagnosticStokesAnalysisEnum(void);
     41int DiagnosticHutterAnalysisEnum(void);
     42int SurfaceSlopeComputeAnalysisEnum(void);
     43int BedSlopeComputeAnalysisEnum(void);
    4144
    4245

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

    r246 r304  
    165165                                        ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
    166166                                        ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
     167                                        ./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
    167168                                        ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
    168169                                        ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
    169170                                        ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
     171                                        ./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
     172                                        ./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
     173                                        ./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
     174                                        ./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
     175                                        ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
     176                                        ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
     177                                        ./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\
     178                                        ./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \
     179                                        ./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\
     180                                        ./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\
     181                                        ./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \
     182                                        ./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\
    170183                                        ./ModelProcessorx/Control/CreateParametersControl.cpp\
    171184                                        ./Dofx/Dofx.h\
    … …  
    187200                                        ./ConfigureObjectsx/ConfigureObjectsx.h\
    188201                                        ./ConfigureObjectsx/ConfigureObjectsx.cpp\
     202                                        ./ComputePressurex/ComputePressurex.h\
     203                                        ./ComputePressurex/ComputePressurex.cpp\
    189204                                        ./BuildNodeSetsx/BuildNodeSetsx.h\
    190205                                        ./BuildNodeSetsx/BuildNodeSetsx.cpp\
    … …  
    387402                                        ./ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp \
    388403                                        ./ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp\
     404                                        ./ModelProcessorx/DiagnosticHoriz/CreateParametersDiagnosticHoriz.cpp\
    389405                                        ./ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp\
    390406                                        ./ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp \
    391407                                        ./ModelProcessorx/DiagnosticVert/CreateLoadsDiagnosticVert.cpp\
     408                                        ./ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp\
     409                                        ./ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \
     410                                        ./ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp\
     411                                        ./ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp\
     412                                        ./ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp \
     413                                        ./ModelProcessorx/DiagnosticStokes/CreateLoadsDiagnosticStokes.cpp\
     414                                        ./ModelProcessorx/SurfaceSlopeCompute/CreateElementsNodesAndMaterialsSurfaceSlopeCompute.cpp\
     415                                        ./ModelProcessorx/SurfaceSlopeCompute/CreateConstraintsSurfaceSlopeCompute.cpp \
     416                                        ./ModelProcessorx/SurfaceSlopeCompute/CreateLoadsSurfaceSlopeCompute.cpp\
     417                                        ./ModelProcessorx/BedSlopeCompute/CreateElementsNodesAndMaterialsBedSlopeCompute.cpp\
     418                                        ./ModelProcessorx/BedSlopeCompute/CreateConstraintsBedSlopeCompute.cpp \
     419                                        ./ModelProcessorx/BedSlopeCompute/CreateLoadsBedSlopeCompute.cpp\
    392420                                        ./ModelProcessorx/Control/CreateParametersControl.cpp\
    393421                                        ./Dofx/Dofx.h\
    … …  
    409437                                        ./ConfigureObjectsx/ConfigureObjectsx.h\
    410438                                        ./ConfigureObjectsx/ConfigureObjectsx.cpp\
     439                                        ./ComputePressurex/ComputePressurex.h\
     440                                        ./ComputePressurex/ComputePressurex.cpp\
    411441                                        ./BuildNodeSetsx/BuildNodeSetsx.h\
    412442                                        ./BuildNodeSetsx/BuildNodeSetsx.cpp\

  • issm/trunk/src/c/NormalizeConstraintsx/NormalizeConstraintsx.cpp

    r1 r304  
    2626        int     null=0;
    2727
    28         /*Recover data: */
    29         msize=nodesets->GetMSize();
    30         nsize=nodesets->GetNSize();
    31         pv_m=nodesets->GetPV_M();
    32         pv_n=nodesets->GetPV_N();
     28        if(nodesets){
     29                /*Recover data: */
     30                msize=nodesets->GetMSize();
     31                nsize=nodesets->GetNSize();
     32                pv_m=nodesets->GetPV_M();
     33                pv_n=nodesets->GetPV_N();
    3334
    3435
    35         if (msize){
     36                if (msize){
    3637
    37                 /*Build row partitioning vector for Rmm and Rmn: */
    38                 row_m=(double*)xmalloc(msize*sizeof(double));
    39                 for(i=0;i<msize;i++)row_m[i]=i+1; //matlab indexing
     38                        /*Build row partitioning vector for Rmm and Rmn: */
     39                        row_m=(double*)xmalloc(msize*sizeof(double));
     40                        for(i=0;i<msize;i++)row_m[i]=i+1; //matlab indexing
    4041
    41                 /*Partition Rmg into Rmm and Rmn: */
    42                 MatPartition(&Rmm,Rmg,row_m,msize,pv_m,msize); //equivalent of Rmm=Rmg(:,mset);
    43                 MatPartition(&Rmn,Rmg,row_m,msize,pv_n,nsize);
     42                        /*Partition Rmg into Rmm and Rmn: */
     43                        MatPartition(&Rmm,Rmg,row_m,msize,pv_m,msize); //equivalent of Rmm=Rmg(:,mset);
     44                        MatPartition(&Rmn,Rmg,row_m,msize,pv_n,nsize);
    4445
    45                 /*Invert Rmm: */
    46                 MatInvert(&InvRmm,Rmm);
     46                        /*Invert Rmm: */
     47                        MatInvert(&InvRmm,Rmm);
    4748
    48                 /*Do Gmn=-(Rmm-1)*Rmn :*/
    49                 MatScale(InvRmm,-1.0);
    50                 MatMatMult(InvRmm,Rmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Gmn);
     49                        /*Do Gmn=-(Rmm-1)*Rmn :*/
     50                        MatScale(InvRmm,-1.0);
     51                        MatMatMult(InvRmm,Rmn,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&Gmn);
    5152
     53                }
     54                else{
     55                        Gmn=NULL;
     56                }
     57
     58                /*Free ressources:*/
     59                xfree((void**)&row_m);
     60                MatFree(&Rmm);
     61                MatFree(&Rmn);
     62                MatFree(&InvRmm);
    5263        }
    53         else{
    54                 Gmn=NULL;
    55         }
    56 
    57         /*Free ressources:*/
    58         xfree((void**)&row_m);
    59         MatFree(&Rmm);
    60         MatFree(&Rmn);
    61         MatFree(&InvRmm);
    6264       
    6365        /*Assign output pointers:*/

  • issm/trunk/src/c/ProcessParamsx/ProcessParamsx.cpp

    r219 r304  
    5151
    5252        /*First serialize partition vector: */
    53         VecToMPISerial(&partition,part);
    54        
     53        if(part)VecToMPISerial(&partition,part);
     54
    5555       
    5656        if (   (analysis_type==ControlAnalysisEnum()) ||
    … …  
    6464                parameters->FindParam((void*)&vz,"vz");
    6565
    66                 u_g=(double*)xcalloc(numberofnodes*3,sizeof(double));
     66                if(numberofnodes){
     67                        u_g=(double*)xcalloc(numberofnodes*3,sizeof(double));
    6768
    68                 if(vx){
    69                         for(i=0;i<numberofnodes;i++){
    70                                 u_g[(int)(3*partition[i]+0)]=vx[i]; 
     69                        if(vx){
     70                                for(i=0;i<numberofnodes;i++){
     71                                        u_g[(int)(3*partition[i]+0)]=vx[i]; 
     72                                }
    7173                        }
    72                 }
    7374
    74                 if(vy){
    75                         for(i=0;i<numberofnodes;i++){
    76                                 u_g[(int)(3*partition[i]+1)]=vy[i]; 
     75                        if(vy){
     76                                for(i=0;i<numberofnodes;i++){
     77                                        u_g[(int)(3*partition[i]+1)]=vy[i]; 
     78                                }
    7779                        }
    78                 }
    79 
    80                 if(vz){
    81                         for(i=0;i<numberofnodes;i++){
    82                                 u_g[(int)(3*partition[i]+2)]=vz[i]; 
     80                        if(vz){
     81                                for(i=0;i<numberofnodes;i++){
     82                                        u_g[(int)(3*partition[i]+2)]=vz[i]; 
     83                                }
    8384                        }
    84                 }
    8585
    8686
    87                 /*Now, create new parameters: */
    88                 count++;
    89                 param= new Param(count,"u_g",DOUBLEVEC);
    90                 param->SetDoubleVec(u_g,3*numberofnodes);
    91                 parameters->AddObject(param);
     87                        /*Now, create new parameters: */
     88                        count++;
     89                        param= new Param(count,"u_g",DOUBLEVEC);
     90                        param->SetDoubleVec(u_g,3*numberofnodes);
     91                        parameters->AddObject(param);
    9292
    93                 /*erase old parameters: */
    94                 param=(Param*)parameters->FindParamObject("vx");
    95                 parameters->DeleteObject((Object*)param);
     93                        /*erase old parameters: */
     94                        param=(Param*)parameters->FindParamObject("vx");
     95                        parameters->DeleteObject((Object*)param);
    9696
    97                 param=(Param*)parameters->FindParamObject("vy");
    98                 parameters->DeleteObject((Object*)param);
     97                        param=(Param*)parameters->FindParamObject("vy");
     98                        parameters->DeleteObject((Object*)param);
    9999
    100                 param=(Param*)parameters->FindParamObject("vz");
    101                 parameters->DeleteObject((Object*)param);
     100                        param=(Param*)parameters->FindParamObject("vz");
     101                        parameters->DeleteObject((Object*)param);
     102                }
    102103
    103104        }

  • issm/trunk/src/c/Reducevectorgtosx/Reducevectorgtosx.cpp

    r1 r304  
    1818        Vec yn=NULL;
    1919
    20         if (nodesets->GetNSize() && nodesets->GetSSize()){
     20        if(nodesets){
    2121
    22                 VecPartition(&yn,yg,nodesets->GetPV_N(),nodesets->GetNSize());
    23        
    24                 VecPartition(&ys,yn,nodesets->GetPV_S(),nodesets->GetSSize());
    25        
     22
     23                if (nodesets->GetNSize() && nodesets->GetSSize()){
     24
     25                        VecPartition(&yn,yg,nodesets->GetPV_N(),nodesets->GetNSize());
     26               
     27                        VecPartition(&ys,yn,nodesets->GetPV_S(),nodesets->GetSSize());
     28               
     29                }
     30
     31                /*Create ys0, full of 0: */
     32                if(ys){
     33                        VecDuplicate(ys,&ys0);
     34                        VecSet(ys0,0.0);
     35                        VecAssemblyBegin(ys0);
     36                        VecAssemblyEnd(ys0);
     37                }
     38
     39                /*Free ressources:*/
     40                VecFree(&yn);
    2641        }
    27 
    28         /*Create ys0, full of 0: */
    29         if(ys){
    30                 VecDuplicate(ys,&ys0);
    31                 VecSet(ys0,0.0);
    32                 VecAssemblyBegin(ys0);
    33                 VecAssemblyEnd(ys0);
    34         }
    35 
    36         /*Free ressources:*/
    37         VecFree(&yn);
    3842       
    3943        /*Assign output pointers:*/

  • issm/trunk/src/c/SpcNodesx/SpcNodesx.cpp

    r1 r304  
    2424        numberofdofs=nodes->NumberOfDofs();
    2525
    26         /*Allocate yg: */
    27         yg=NewVec(numberofdofs);
     26        if(numberofdofs){
     27                /*Allocate yg: */
     28                yg=NewVec(numberofdofs);
    2829
    29         /*Now, go through constraints, and update the nodes and the constraint vector at the same time: */
    30         constraints->SetupSpcs(nodes,yg);
     30                /*Now, go through constraints, and update the nodes and the constraint vector at the same time: */
     31                constraints->SetupSpcs(nodes,yg);
     32        }
    3133
    3234        /*Assign output pointers: */

  • issm/trunk/src/c/io/FetchNodeSets.cpp

    r1 r304  
    3232        int ssize;
    3333
    34         FetchData((void**)&pv_m,NULL,NULL,mxGetField(dataref,0,"pv_m"),"Vector","Vec");
    35         FetchData((void**)&pv_n,NULL,NULL,mxGetField(dataref,0,"pv_n"),"Vector","Vec");
    36         FetchData((void**)&pv_f,NULL,NULL,mxGetField(dataref,0,"pv_f"),"Vector","Vec");
    37         FetchData((void**)&pv_s,NULL,NULL,mxGetField(dataref,0,"pv_s"),"Vector","Vec");
    38        
    39         gsize=(int)mxGetScalar(mxGetField(dataref,0,"gsize"));
    40         msize=(int)mxGetScalar(mxGetField(dataref,0,"msize"));
    41         nsize=(int)mxGetScalar(mxGetField(dataref,0,"nsize"));
    42         fsize=(int)mxGetScalar(mxGetField(dataref,0,"fsize"));
    43         ssize=(int)mxGetScalar(mxGetField(dataref,0,"ssize"));
     34        if(mxIsEmpty(dataref)){
     35                nodesets=NULL;
     36        }
     37        else{
    4438
    45         nodesets=new NodeSets( pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
     39                FetchData((void**)&pv_m,NULL,NULL,mxGetField(dataref,0,"pv_m"),"Vector","Vec");
     40                FetchData((void**)&pv_n,NULL,NULL,mxGetField(dataref,0,"pv_n"),"Vector","Vec");
     41                FetchData((void**)&pv_f,NULL,NULL,mxGetField(dataref,0,"pv_f"),"Vector","Vec");
     42                FetchData((void**)&pv_s,NULL,NULL,mxGetField(dataref,0,"pv_s"),"Vector","Vec");
     43               
     44                gsize=(int)mxGetScalar(mxGetField(dataref,0,"gsize"));
     45                msize=(int)mxGetScalar(mxGetField(dataref,0,"msize"));
     46                nsize=(int)mxGetScalar(mxGetField(dataref,0,"nsize"));
     47                fsize=(int)mxGetScalar(mxGetField(dataref,0,"fsize"));
     48                ssize=(int)mxGetScalar(mxGetField(dataref,0,"ssize"));
     49
     50                nodesets=new NodeSets( pv_m,pv_n,pv_f,pv_s,gsize,msize,nsize,fsize,ssize);
     51        }
    4652
    4753        /*Assign output pointers:*/

  • issm/trunk/src/c/io/WriteNodeSets.cpp

    r1 r304  
    4141        mxArray*    field=NULL;
    4242
    43         /*Recover data from the nodesets class: */
    44         gsize=nodesets->GetGSize();
    45         msize=nodesets->GetMSize();
    46         nsize=nodesets->GetNSize();
    47         fsize=nodesets->GetFSize();
    48         ssize=nodesets->GetSSize();
     43        if(nodesets){
     44                /*Recover data from the nodesets class: */
     45                gsize=nodesets->GetGSize();
     46                msize=nodesets->GetMSize();
     47                nsize=nodesets->GetNSize();
     48                fsize=nodesets->GetFSize();
     49                ssize=nodesets->GetSSize();
    4950
    50         if(msize){
    51                 pv_m=(double*)xmalloc(msize*sizeof(double));
    52                 memcpy(pv_m,nodesets->GetPV_M(),msize*sizeof(double));
    53         }
    54         if(nsize){
    55                 pv_n=(double*)xmalloc(nsize*sizeof(double));
    56                 memcpy(pv_n,nodesets->GetPV_N(),nsize*sizeof(double));
    57         }
    58         if(fsize){
    59                 pv_f=(double*)xmalloc(fsize*sizeof(double));
    60                 memcpy(pv_f,nodesets->GetPV_F(),fsize*sizeof(double));
    61         }
    62         if(ssize){
    63                 pv_s=(double*)xmalloc(ssize*sizeof(double));
    64                 memcpy(pv_s,nodesets->GetPV_S(),ssize*sizeof(double));
    65         }
     51                if(msize){
     52                        pv_m=(double*)xmalloc(msize*sizeof(double));
     53                        memcpy(pv_m,nodesets->GetPV_M(),msize*sizeof(double));
     54                }
     55                if(nsize){
     56                        pv_n=(double*)xmalloc(nsize*sizeof(double));
     57                        memcpy(pv_n,nodesets->GetPV_N(),nsize*sizeof(double));
     58                }
     59                if(fsize){
     60                        pv_f=(double*)xmalloc(fsize*sizeof(double));
     61                        memcpy(pv_f,nodesets->GetPV_F(),fsize*sizeof(double));
     62                }
     63                if(ssize){
     64                        pv_s=(double*)xmalloc(ssize*sizeof(double));
     65                        memcpy(pv_s,nodesets->GetPV_S(),ssize*sizeof(double));
     66                }
    6667
    67         /*Build structure in matlab workspace with all these fields: */
    68         fnames[0] = "gsize";
    69         fnames[1] = "msize";
    70         fnames[2] = "nsize";
    71         fnames[3] = "fsize";
    72         fnames[4] = "ssize";
    73         fnames[5] = "pv_m";
    74         fnames[6] = "pv_n";
    75         fnames[7] = "pv_f";
    76         fnames[8] = "pv_s";
     68                /*Build structure in matlab workspace with all these fields: */
     69                fnames[0] = "gsize";
     70                fnames[1] = "msize";
     71                fnames[2] = "nsize";
     72                fnames[3] = "fsize";
     73                fnames[4] = "ssize";
     74                fnames[5] = "pv_m";
     75                fnames[6] = "pv_n";
     76                fnames[7] = "pv_f";
     77                fnames[8] = "pv_s";
    7778
    78         dataref=mxCreateStructArray( ndim,onebyone,nfields,fnames);
    79        
    80         mxSetField( dataref, 0, "gsize",mxCreateDoubleScalar((double)gsize));
    81         mxSetField( dataref, 0, "msize",mxCreateDoubleScalar((double)msize));
    82         mxSetField( dataref, 0, "nsize",mxCreateDoubleScalar((double)nsize));
    83         mxSetField( dataref, 0, "fsize",mxCreateDoubleScalar((double)fsize));
    84         mxSetField( dataref, 0, "ssize",mxCreateDoubleScalar((double)ssize));
     79                dataref=mxCreateStructArray( ndim,onebyone,nfields,fnames);
     80               
     81                mxSetField( dataref, 0, "gsize",mxCreateDoubleScalar((double)gsize));
     82                mxSetField( dataref, 0, "msize",mxCreateDoubleScalar((double)msize));
     83                mxSetField( dataref, 0, "nsize",mxCreateDoubleScalar((double)nsize));
     84                mxSetField( dataref, 0, "fsize",mxCreateDoubleScalar((double)fsize));
     85                mxSetField( dataref, 0, "ssize",mxCreateDoubleScalar((double)ssize));
    8586
    86         if(msize){
    87                 field = mxCreateDoubleMatrix(0,0,mxREAL);
    88                 mxSetM(field,msize); mxSetN(field,1); mxSetPr(field,pv_m);
     87                if(msize){
     88                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     89                        mxSetM(field,msize); mxSetN(field,1); mxSetPr(field,pv_m);
     90                }
     91                else{
     92                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     93                }
     94                mxSetField( dataref, 0, "pv_m",field);
     95
     96                if(nsize){
     97                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     98                        mxSetM(field,nsize); mxSetN(field,1); mxSetPr(field,pv_n);
     99                }
     100                else{
     101                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     102                }
     103                mxSetField( dataref, 0, "pv_n",field);
     104
     105               
     106                if(fsize){
     107                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     108                        mxSetM(field,fsize); mxSetN(field,1); mxSetPr(field,pv_f);
     109                }
     110                else{
     111                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     112                }
     113                mxSetField( dataref, 0, "pv_f",field);
     114               
     115                if(ssize){
     116                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     117                        mxSetM(field,ssize); mxSetN(field,1); mxSetPr(field,pv_s);
     118                }
     119                else{
     120                        field = mxCreateDoubleMatrix(0,0,mxREAL);
     121                }
     122                mxSetField( dataref, 0, "pv_s",field);
    89123        }
    90124        else{
    91                 field = mxCreateDoubleMatrix(0,0,mxREAL);
     125                dataref=mxCreateDoubleMatrix(0,0,mxREAL);
    92126        }
    93         mxSetField( dataref, 0, "pv_m",field);
    94 
    95         if(nsize){
    96                 field = mxCreateDoubleMatrix(0,0,mxREAL);
    97                 mxSetM(field,nsize); mxSetN(field,1); mxSetPr(field,pv_n);
    98         }
    99         else{
    100                 field = mxCreateDoubleMatrix(0,0,mxREAL);
    101         }
    102         mxSetField( dataref, 0, "pv_n",field);
    103 
    104        
    105         if(fsize){
    106                 field = mxCreateDoubleMatrix(0,0,mxREAL);
    107                 mxSetM(field,fsize); mxSetN(field,1); mxSetPr(field,pv_f);
    108         }
    109         else{
    110                 field = mxCreateDoubleMatrix(0,0,mxREAL);
    111         }
    112         mxSetField( dataref, 0, "pv_f",field);
    113        
    114         if(ssize){
    115                 field = mxCreateDoubleMatrix(0,0,mxREAL);
    116                 mxSetM(field,ssize); mxSetN(field,1); mxSetPr(field,pv_s);
    117         }
    118         else{
    119                 field = mxCreateDoubleMatrix(0,0,mxREAL);
    120         }
    121         mxSetField( dataref, 0, "pv_s",field);
    122127
    123128        /*Assign output pointers:*/

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

    r132 r304  
    4646#include "./VelocityExtrudex/VelocityExtrudex.h"
    4747#include "./GriddataMeshToGridx/GriddataMeshToGridx.h"
     48#include "./ComputePressurex/ComputePressurex.h"
    4849
    4950

  • issm/trunk/src/c/shared/Numerics/DistributeNumDofs.cpp

    r128 r304  
    2222                numdofs=1;
    2323        }
     24        else if (analysis_type==DiagnosticStokesAnalysisEnum()){
     25                numdofs=4;
     26        }
     27        else if (analysis_type==DiagnosticHutterAnalysisEnum()){
     28                numdofs=2;
     29        }
     30        else if (analysis_type==SurfaceSlopeComputeAnalysisEnum()){
     31                numdofs=1;
     32        }
     33        else if (analysis_type==BedSlopeComputeAnalysisEnum()){
     34                numdofs=1;
     35        }
    2436        else if (analysis_type==ControlAnalysisEnum()){
    2537                numdofs=2;

  • issm/trunk/src/c/toolkits/petsc/patches/PetscVectorToMatlabVector.cpp

    r1 r304  
    3939        VecGetSize(vector,&rows);
    4040
    41         idxm=(int*)xmalloc(rows*sizeof(int));
    42         values=(double*)xmalloc(rows*sizeof(double));
    43         for(i=0;i<rows;i++)idxm[i]=i;
    44          
    45         VecGetValues(vector,rows,idxm,values);
    46 
     41        if(rows){
     42                idxm=(int*)xmalloc(rows*sizeof(int));
     43                values=(double*)xmalloc(rows*sizeof(double));
     44                for(i=0;i<rows;i++)idxm[i]=i;
     45                 
     46                VecGetValues(vector,rows,idxm,values);
     47        }
    4748
    4849        /*Using values, build a matlab vector: */

  • issm/trunk/src/m/classes/public/setelementstype.m

    r202 r304  
    121121
    122122%figure out solution types
    123 md.ishutter=any(md.elements_type(:,1)==hutterenum());
    124 md.ismacayealpattyn=any(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum() );
    125 md.isstokes=any(md.elements_type(:,2));
     123md.ishutter=double(any(md.elements_type(:,1)==hutterenum()));
     124md.ismacayealpattyn=double(any(md.elements_type(:,1)==macayealenum() | md.elements_type(:,1)==pattynenum() ));
     125md.isstokes=double(any(md.elements_type(:,2)));
    126126
    127127end

  • issm/trunk/src/m/classes/public/solve.m

    r152 r304  
    5454end
    5555
    56 disp(sprintf('\nlaunching solution sequence\n'));
     56displaystring(md.debug,'%s','\nlaunching solution sequence\n');
    5757
    5858%If running in parallel, we have a different way of launching the solution
    … …  
    9797
    9898%Check result is consistent
    99 disp(sprintf('%s\n','checking result consistency'));
     99displaystring(md.debug,'%s\n','checking result consistency');
    100100if ~isresultconsistent(md,solutiontype),
    101101        disp('!! results not consistent correct the model !!') %it would be very cruel to put an error, it would kill the computed results (even if not consistent...)

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

    r132 r304  
    66else
    77bin_PROGRAMS =  BuildNodeSets\
     8                                ComputePressure\
    89                                ConfigureObjects \
    910                                ControlOptimization\
    … …  
    6869BuildNodeSets_SOURCES = BuildNodeSets/BuildNodeSets.cpp\
    6970                          BuildNodeSets/BuildNodeSets.h
     71
     72ComputePressure_SOURCES = ComputePressure/ComputePressure.cpp\
     73                          ComputePressure/ComputePressure.h
    7074
    7175ConfigureObjects_SOURCES = ConfigureObjects/ConfigureObjects.cpp\

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

    r200 r304  
    2525
    2626        /*Shit partition: */
    27         VecShift(partition,-1.0);
     27        if(partition)VecShift(partition,-1.0);
    2828
    2929        /*!Generate internal degree of freedom numbers: */

  • issm/trunk/todo

    r226 r304  
    66Update control to get new inputs, as well as parallel diagnostic
    77
     8Finish diagnostic_core.m
Note: See TracChangeset for help on using the changeset viewer.