Changeset 586


Ignore:
Timestamp:
05/26/09 10:54:45 (15 years ago)
Author:
Eric.Larour
Message:

Added dakota parallel c code capability + Prognostic capability

Location:
issm/trunk/src
Files:
31 added
42 edited

Legend:

Unmodified
Added
Removed

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

    r503 r586  
    12291229
    12301230}
    1231 
     1231void  DataSet::ThicknessExtrude(Vec tg,double* tg_serial){
     1232
     1233        vector<Object*>::iterator object;
     1234        Penta* penta=NULL;
     1235
     1236        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1237               
     1238                if((*object)->Enum()==PentaEnum()){
     1239
     1240                        penta=(Penta*)(*object);
     1241                        penta->ThicknessExtrude(tg,tg_serial);
     1242
     1243                }
     1244        }
     1245
     1246}
     1247void  DataSet::VelocityExtrudeAllElements(Vec ug,double* ug_serial){
     1248
     1249        vector<Object*>::iterator object;
     1250        Penta* penta=NULL;
     1251
     1252        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1253               
     1254                if((*object)->Enum()==PentaEnum()){
     1255
     1256                        penta=(Penta*)(*object);
     1257                        penta->VelocityExtrudeAllElements(ug,ug_serial);
     1258
     1259                }
     1260        }
     1261
     1262}
     1263void  DataSet::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){
     1264
     1265       
     1266        vector<Object*>::iterator object;
     1267        Penta* penta=NULL;
     1268
     1269        for ( object=objects.begin() ; object < objects.end(); object++ ){
     1270               
     1271                if((*object)->Enum()==PentaEnum()){
     1272
     1273                        penta=(Penta*)(*object);
     1274                        penta->VelocityDepthAverageAtBase(ug,ug_serial);
     1275
     1276                }
     1277        }
     1278
     1279}
    12321280void  DataSet::SlopeExtrude(Vec sg,double* sg_serial){
    12331281

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

    r465 r586  
    7575                void  Misfit(double* pJ, double* u_g,double* u_g_obs,void* inputs,int analysis_type,int sub_analysis_type);
    7676                void  VelocityExtrude(Vec ug,double* ug_serial);
     77                void  ThicknessExtrude(Vec ug,double* ug_serial);
     78                void  VelocityExtrudeAllElements(Vec ug,double* ug_serial);
     79                void  VelocityDepthAverageAtBase(Vec ug,double* ug_serial);
    7780                void  SlopeExtrude(Vec sg,double* sg_serial);
    7881                int   DeleteObject(Object* object);

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

    r483 r586  
    1 INCLUDES = @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@  @METISINCL@ @PLAPACKINCL@  @BLASLAPACKINCL@ @MUMPSINCL@  @TRIANGLEINCL@
     1INCLUDES = @DAKOTAINCL@ @PETSCINCL@ @SLEPCINCL@ @MPIINCL@ @MATLABINCL@  @METISINCL@ @PLAPACKINCL@  @BLASLAPACKINCL@ @MUMPSINCL@  @TRIANGLEINCL@
    22
    33#Compile serial library, and then try and compile parallel library
    … …  
    126126                                        ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\
    127127                                        ./toolkits/petsc/patches/NewMat.cpp\
     128                                        ./toolkits/petsc/patches/SerialToVec.cpp\
    128129                                        ./toolkits/petsc/patches/VecFree.cpp\
    129130                                        ./toolkits/petsc/patches/VecDuplicatePatch.cpp\
    … …  
    196197                                        ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\
    197198                                        ./ModelProcessorx/Melting/CreateParametersMelting.cpp\
     199                                        ./ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp\
     200                                        ./ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp\
     201                                        ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\
     202                                        ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\
    198203                                        ./Dofx/Dofx.h\
    199204                                        ./Dofx/Dofx.cpp\
    … …  
    252257                                        ./ProcessParamsx/ProcessParamsx.cpp\
    253258                                        ./ProcessParamsx/ProcessParamsx.h\
     259                                        ./ThicknessExtrudex/ThicknessExtrudex.cpp\
     260                                        ./ThicknessExtrudex/ThicknessExtrudex.h\
    254261                                        ./VelocityExtrudex/VelocityExtrudex.cpp\
    255262                                        ./VelocityExtrudex/VelocityExtrudex.h\
     263                                        ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\
     264                                        ./VelocityDepthAveragex/VelocityDepthAveragex.h\
    256265                                        ./SlopeExtrudex/SlopeExtrudex.cpp\
    257266                                        ./SlopeExtrudex/SlopeExtrudex.h
    … …  
    288297                                        ./objects/Node.h\
    289298                                        ./objects/Node.cpp\
     299                                        ./objects/DakotaPlugin.h\
     300                                        ./objects/DakotaPlugin.cpp\
    290301                                        ./objects/Tria.h\
    291302                                        ./objects/Tria.cpp\
    … …  
    380391                                        ./toolkits/petsc/patches/PetscOptionsInsertMultipleString.cpp\
    381392                                        ./toolkits/petsc/patches/NewMat.cpp\
     393                                        ./toolkits/petsc/patches/SerialToVec.cpp\
    382394                                        ./toolkits/petsc/patches/VecFree.cpp\
    383395                                        ./toolkits/petsc/patches/VecDuplicatePatch.cpp\
    … …  
    450462                                        ./ModelProcessorx/Melting/CreateLoadsMelting.cpp\
    451463                                        ./ModelProcessorx/Melting/CreateParametersMelting.cpp\
     464                                        ./ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp\
     465                                        ./ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp\
     466                                        ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\
     467                                        ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\
    452468                                        ./Dofx/Dofx.h\
    453469                                        ./Dofx/Dofx.cpp\
    … …  
    502518                                        ./Solverx/Solverx.cpp\
    503519                                        ./Solverx/Solverx.h\
     520                                        ./ThicknessExtrudex/ThicknessExtrudex.cpp\
     521                                        ./ThicknessExtrudex/ThicknessExtrudex.h\
    504522                                        ./Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\
    505523                                        ./Mergesolutionfromftogx/Mergesolutionfromftogx.h\
    … …  
    508526                                        ./VelocityExtrudex/VelocityExtrudex.cpp\
    509527                                        ./VelocityExtrudex/VelocityExtrudex.h\
     528                                        ./VelocityDepthAveragex/VelocityDepthAveragex.cpp\
     529                                        ./VelocityDepthAveragex/VelocityDepthAveragex.h\
    510530                                        ./SlopeExtrudex/SlopeExtrudex.cpp\
    511531                                        ./SlopeExtrudex/SlopeExtrudex.h\
    … …  
    520540                                        ./parallel/OutputControl.cpp\
    521541                                        ./parallel/OutputThermal.cpp\
     542                                        ./parallel/OutputPrognostic.cpp\
    522543                                        ./parallel/objectivefunctionC.cpp\
    523                                         ./parallel/GradJCompute.cpp
     544                                        ./parallel/GradJCompute.cpp\
     545                                        ./parallel/qmu.cpp\
     546                                        ./parallel/SpawnCore.cpp\
     547                                        ./parallel/DakotaResponses.cpp
    524548
    525549libpISSM_a_CXXFLAGS = -fPIC -D_PARALLEL_   -D_C_
    … …  
    528552bin_PROGRAMS =
    529553else
    530 bin_PROGRAMS = diagnostic.exe control.exe thermal.exe
     554bin_PROGRAMS = diagnostic.exe
     555#control.exe thermal.exe prognostic.exe
    531556endif
    532557
    533 LDADD = ./libpISSM.a $(TRIANGLELIB) $(METISLIB) $(PETSCLIB) $(SLEPCLIB) $(MUMPSLIB) $(PLAPACKLIB)  $(MPILIB) $(X_LIBS) -lX11 $(BLASLAPACKLIB) $(SCALAPACKLIB) $(BLACSLIB) $(FLIBS)
     558LDADD = ./libpISSM.a $(TRIANGLELIB) $(METISLIB) $(PETSCLIB) $(DAKOTALIB) $(SLEPCLIB) $(MUMPSLIB) $(PLAPACKLIB)  $(MPILIB) $(X_LIBS) -lX11 $(BLASLAPACKLIB) $(SCALAPACKLIB) $(BLACSLIB) $(FLIBS)
    534559
    535560diagnostic_exe_SOURCES = parallel/diagnostic.cpp
    … …  
    541566thermal_exe_SOURCES = parallel/thermal.cpp
    542567thermal_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     568
     569prognostic_exe_SOURCES = parallel/prognostic.cpp
     570prognostic_exe_CXXFLAGS= -fPIC -D_PARALLEL_
     571

  • issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp

    r576 r586  
    8787        }
    8888        else if (strcmp(model->analysis_type,"prognostic")==0){
    89                        
    90                 //CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, model,model_handle);
    91                 //CreateConstraintsPrognostic(pconstraints,model,model_handle);
    92                 //CreateLoadsPrognostic(ploads,model,model_handle);
     89
     90                CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pmaterials, model,model_handle);
     91                CreateConstraintsPrognostic(pconstraints,model,model_handle);
     92                CreateLoadsPrognostic(ploads,model,model_handle);
     93                CreateParametersPrognostic(pparameters,model,model_handle);
     94                                       
    9395        }
    9496        else{

  • issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp

    r465 r586  
    2424        int      numberofdofspernode;
    2525        int      dim;
     26        double*  epart=NULL;
     27        double*  part=NULL;
     28
    2629
    2730        /*Initialize dataset: */
    … …  
    4245        parameters->AddObject(param);
    4346
     47        //qmu analysis?
     48        count++;
     49        param= new Param(count,"qmu_analysis",INTEGER);
     50        param->SetInteger(model->qmu_analysis);
     51        parameters->AddObject(param);
     52
     53        count++;
     54        param= new Param(count,"qmu_npart",INTEGER);
     55        param->SetInteger(model->qmu_npart);
     56        parameters->AddObject(param);
     57
    4458        //dimension 2d or 3d:
    4559        if (strcmp(model->meshtype,"2d")==0)dim=2;
    … …  
    201215        parameters->AddObject(param);
    202216
     217        /*Deal with responses and partition for qmu modeling: */
     218        if(model->qmu_analysis){
     219                char** descriptors=NULL;
     220                char*  descriptor=NULL;
     221                char* tag=NULL;
     222
     223                descriptors=(char**)xmalloc(model->numberofresponses*sizeof(char*));
     224                tag=(char*)xmalloc((strlen("descriptori")+1)*sizeof(char));
     225
     226                /*Fetch descriptors: */
     227                for(i=0;i<model->numberofresponses;i++){
     228                        sprintf(tag,"%s%i","descriptor",i);
     229                        ModelFetchData((void**)&descriptor,NULL,NULL,model_handle,tag,"String",NULL);
     230                        descriptors[i]=descriptor;
     231                }
     232
     233                /*Ok, we have all the descriptors. Build a parameter with it: */
     234                count++;
     235                param= new Param(count,"descriptors",STRINGARRAY);
     236                param->SetStringArray(descriptors,model->numberofresponses);
     237                parameters->AddObject(param);
     238
     239                /*Free data: */
     240                xfree((void**)&tag);
     241                for(i=0;i<model->numberofresponses;i++){
     242                        char* descriptor=descriptors[i];
     243                        xfree((void**)&descriptor);
     244                }
     245                xfree((void**)&descriptors);
     246
     247                #ifdef _PARALLEL_
     248                /*partition grids in model->qmu_npart parts: */
     249       
     250                if(strcmp(model->meshtype,"2d")==0){
     251                        ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat");
     252                }
     253                else{
     254                        ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat");
     255                }
     256
     257                MeshPartitionx(&epart, &part,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,model->qmu_npart);
     258
     259                count++;
     260                param= new Param(count,"qmu_part",DOUBLEVEC);
     261                param->SetDoubleVec(part,model->numberofnodes);
     262                parameters->AddObject(param);
     263
     264                /*Free elements and elements2d: */
     265                xfree((void**)&model->elements);
     266                xfree((void**)&model->elements2d);
     267                xfree((void**)&epart);
     268                xfree((void**)&part);
     269
     270                #endif
     271        }
    203272
    204273       

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp

    r300 r586  
    8484
    8585       
    86         cleanup_and_return:
    8786        /*Free data: */
    8887        xfree((void**)&gridondirichlet_diag);
    8988        xfree((void**)&dirichletvalues_diag);
    9089       
     90        cleanup_and_return:
    9191        /*Assign output pointer: */
    9292        *pconstraints=constraints;

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp

    r483 r586  
    6969        double tria_meanvel;/*!scaling ratio for velocities*/
    7070        double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
    71         double tria_viscosity_overshoot;
    72 
     71        double tria_viscosity_overshoot;
     72        int    tria_artdiff;
     73       
    7374        /*matice constructor input: */
    7475        int    matice_mid;
    … …  
    294295
    295296                        /*Create tria element using its constructor:*/
    296                         tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
     297                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff);
    297298
    298299                        /*Add tria element to elements dataset: */
    … …  
    653654        xfree((void**)&model->uppernodes);
    654655               
    655         cleanup_and_return:
    656 
    657         /*Assign output pointer: */
    658         *pelements=elements;
    659         *pnodes=nodes;
    660         *pmaterials=materials;
    661656
    662657        /*Keep partitioning information into model*/
    … …  
    672667        #endif
    673668
     669        cleanup_and_return:
     670
     671        /*Assign output pointer: */
     672        *pelements=elements;
     673        *pnodes=nodes;
     674        *pmaterials=materials;
     675
    674676}

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp

    r300 r586  
    265265        loads->Presort();
    266266
    267         cleanup_and_return:
    268267
    269268        /*Free ressources:*/
    … …  
    276275        xfree((void**)&riftsfill);
    277276        xfree((void**)&riftsfriction);
    278 
     277       
     278        cleanup_and_return:
    279279
    280280        /*Assign output pointer: */

  • issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp

    r465 r586  
    357357        xfree((void**)&model->uppernodes);
    358358               
     359
     360        /*Keep partitioning information into model*/
     361        xfree((void**)&epart);
     362        model->npart=npart;
     363
    359364        cleanup_and_return:
    360365
    … …  
    364369        *pmaterials=materials;
    365370
    366         /*Keep partitioning information into model*/
    367         xfree((void**)&epart);
    368         model->npart=npart;
    369371}

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateConstraintsDiagnosticStokes.cpp

    r394 r586  
    8888        constraints->Presort();
    8989
    90         cleanup_and_return:
    9190        /*Free data: */
    9291        xfree((void**)&gridonstokes);
     92
     93        cleanup_and_return:
    9394
    9495        /*Assign output pointer: */

  • issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp

    r465 r586  
    517517        xfree((void**)&model->borderstokes);
    518518
    519         cleanup_and_return:
    520 
    521         /*Assign output pointer: */
    522         *pelements=elements;
    523         *pnodes=nodes;
    524         *pmaterials=materials;
    525519
    526520        /*Keep partitioning information into model*/
    … …  
    536530        #endif
    537531
     532        cleanup_and_return:
     533
     534        /*Assign output pointer: */
     535        *pelements=elements;
     536        *pnodes=nodes;
     537        *pmaterials=materials;
     538
    538539}

  • issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp

    r465 r586  
    146146        sub_analysis_type=AnalysisTypeAsEnum(model->sub_analysis_type);
    147147
    148         /*Width of elements: */
    149         if(strcmp(model->meshtype,"2d")==0){
    150                 throw ErrorException(__FUNCT__," 2d mesh not supported for vertical velocity");
    151         }
    152 
    153 
    154148        #ifdef _PARALLEL_
    155149        /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
    … …  
    420414        xfree((void**)&model->uppernodes);
    421415               
     416
     417        /*Keep partitioning information into model*/
     418        model->epart=epart;
     419        model->my_grids=my_grids;
     420        model->my_bordergrids=my_bordergrids;
     421
     422        /*Free ressources:*/
     423        #ifdef _PARALLEL_
     424        xfree((void**)&all_numgrids);
     425        xfree((void**)&npart);
     426        VecFree(&gridborder);
     427        #endif
     428
    422429        cleanup_and_return:
    423430
    … …  
    427434        *pmaterials=materials;
    428435
    429         /*Keep partitioning information into model*/
    430         model->epart=epart;
    431         model->my_grids=my_grids;
    432         model->my_bordergrids=my_bordergrids;
    433 
    434         /*Free ressources:*/
    435         #ifdef _PARALLEL_
    436         xfree((void**)&all_numgrids);
    437         xfree((void**)&npart);
    438         VecFree(&gridborder);
    439         #endif
    440 
    441436}

  • issm/trunk/src/c/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp

    r483 r586  
    476476        xfree((void**)&model->uppernodes);
    477477               
     478
     479        /*Keep partitioning information into model*/
     480        model->epart=epart;
     481        model->my_grids=my_grids;
     482        model->my_bordergrids=my_bordergrids;
     483
     484        /*Free ressources:*/
     485        #ifdef _PARALLEL_
     486        xfree((void**)&all_numgrids);
     487        xfree((void**)&npart);
     488        VecFree(&gridborder);
     489        #endif
     490
    478491        cleanup_and_return:
    479492
    … …  
    483496        *pmaterials=materials;
    484497
    485         /*Keep partitioning information into model*/
    486         model->epart=epart;
    487         model->my_grids=my_grids;
    488         model->my_bordergrids=my_bordergrids;
    489 
    490         /*Free ressources:*/
    491         #ifdef _PARALLEL_
    492         xfree((void**)&all_numgrids);
    493         xfree((void**)&npart);
    494         VecFree(&gridborder);
    495         #endif
    496 
    497498}

  • issm/trunk/src/c/ModelProcessorx/Model.cpp

    r575 r586  
    3232
    3333        /*!initialize all pointers to 0: */
    34 
    35 
    3634        model->repository=NULL;
    3735        model->meshtype=NULL;
    3836        model->analysis_type=NULL;
    3937        model->sub_analysis_type=NULL;
     38        model->qmu_analysis=0;
     39        model->solverstring=NULL;
     40        model->numberofresponses=0;
     41        model->qmu_npart=0;
    4042        model->numberofelements=0;
    4143        model->numberofnodes=0;
    … …  
    7577        model->drag=NULL;
    7678        model->p=NULL;
     79        model->q=NULL;
    7780       
    7881       
    … …  
    100103        model->rho_water=0;
    101104        model->rho_ice=0;
    102                 model->g=0;
    103         model->n=0;
     105        model->g=0;
     106        model->n=NULL;
    104107        model->B=NULL;
    105108
    … …  
    165168        model->ismacayealpattyn=0;
    166169        model->isstokes=0;
     170
     171
     172        model->epart=NULL;
     173        model->npart=NULL;
     174        model->my_grids=NULL;
     175        model->my_bordergrids=NULL;
    167176
    168177        return model;
    … …  
    295304
    296305        /*In ModelInit, we get all the data that is not difficult to get, and that is small: */
     306
    297307        ModelFetchData((void**)&model->analysis_type,NULL,NULL,model_handle,"analysis_type","String",NULL);
    298308        ModelFetchData((void**)&model->sub_analysis_type,NULL,NULL,model_handle,"sub_analysis_type","String",NULL);
     309        ModelFetchData((void**)&model->qmu_analysis,NULL,NULL,model_handle,"qmu_analysis","Integer",NULL);
    299310
    300311        ModelFetchData((void**)&model->meshtype,NULL,NULL,model_handle,"type","String",NULL);
    … …  
    302313        ModelFetchData((void**)&model->numberofnodes,NULL,NULL,model_handle,"numberofgrids","Integer",NULL);
    303314        ModelFetchData((void**)&model->numberofelements,NULL,NULL,model_handle,"numberofelements","Integer",NULL);
    304 
    305315        /*!In case we are running 3d, we are going to need the collapsed and non-collapsed 2d meshes, from which the 3d mesh was extruded: */
    306316        if (strcmp(model->meshtype,"3d")==0){
    … …  
    311321                ModelFetchData((void**)&model->numlayers,NULL,NULL,model_handle,"numlayers","Integer",NULL);
    312322        }
     323
    313324
    314325        /*elements type: */
    … …  
    367378        ModelFetchData((void**)&model->numrifts,NULL,NULL,model_handle,"numrifts","Integer",NULL);
    368379
     380        /*qmu: */
     381        if(model->qmu_analysis){
     382                ModelFetchData((void**)&model->numberofresponses,NULL,NULL,model_handle,"numberofresponses","Integer",NULL);
     383                ModelFetchData((void**)&model->qmu_npart,NULL,NULL,model_handle,"npart","Integer",NULL);
     384        }
     385
    369386        /*Assign output pointers: */
    370387        *pmodel=model;

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

    r575 r586  
    1717        char*   analysis_type;
    1818        char*   sub_analysis_type;
     19        int     qmu_analysis;
    1920        char*   solverstring;
    2021
    … …  
    5253        double*  vx_obs;
    5354        double*  vy_obs;
     55
     56        /*qmu: */
     57        int      numberofresponses;
     58        int      qmu_npart;
    5459
    5560        /*geometry: */
    … …  
    221226        void    CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
    222227        void    CreateConstraintsMelting(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
    223         void  CreateLoadsMelting(DataSet** ploads, Model* model, ConstDataHandle model_handle);
    224         void  CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
     228        void    CreateLoadsMelting(DataSet** ploads, Model* model, ConstDataHandle model_handle);
     229        void    CreateParametersMelting(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
     230
     231        /*prognostic:*/
     232        void    CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle);
     233        void    CreateConstraintsPrognostic(DataSet** pconstraints,Model* model,ConstDataHandle model_handle);
     234        void    CreateLoadsPrognostic(DataSet** ploads, Model* model, ConstDataHandle model_handle);
     235        void    CreateParametersPrognostic(DataSet** pparameters,Model* model,ConstDataHandle model_handle);
    225236
    226237

  • issm/trunk/src/c/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp

    r483 r586  
    6666        double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/
    6767        double tria_viscosity_overshoot;
     68        int    tria_artdiff;
    6869
    6970        /*penta constructor input: */
    … …  
    227228
    228229                        /*Create tria element using its constructor:*/
    229                         tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot);
     230                        tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff);
    230231
    231232                        /*Add tria element to elements dataset: */
    … …  
    479480        xfree((void**)&model->uppernodes);
    480481               
    481         cleanup_and_return:
    482 
    483         /*Assign output pointer: */
    484         *pelements=elements;
    485         *pnodes=nodes;
    486         *pmaterials=materials;
    487482
    488483        /*Keep partitioning information into model*/
    … …  
    497492        VecFree(&gridborder);
    498493        #endif
     494
     495        cleanup_and_return:
     496
     497        /*Assign output pointer: */
     498        *pelements=elements;
     499        *pnodes=nodes;
     500        *pmaterials=materials;
     501
    499502}

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

    r543 r586  
    4747        }
    4848
     49        #ifdef _DEBUG_
     50        MatNorm(Kgg,NORM_INFINITY,&kmax2);
     51         _printf_("   K_gg infinity norm after penalties: %g\n",kmax2);
     52        #endif
     53
    4954        /*Assign output pointers:*/
    5055        if(pkmax)*pkmax=kmax;

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

    r552 r586  
    2424        double* vy=NULL;
    2525        double* vz=NULL;
    26        
     26        double* pressure=NULL;
     27               
    2728        double* u_g=NULL;
    28 
    29         /*thermal*/
    30         double* pressure=NULL;
    31         double* temperature=NULL;
    32         double* melting=NULL;
    33 
    34         double* p_g=NULL;
    35         double* t_g=NULL;
    36         double* m_g=NULL;
    3729       
    3830        /*control: */
    … …  
    4638
    4739        double* u_g_obs=NULL;
     40        double* p_g=NULL;
     41
     42        /*prognostic: */
     43        double* a_g=NULL;
     44        double* m_g=NULL;
     45        double* h_g=NULL;
     46        double* accumulation=NULL;
     47        double* thickness=NULL;
     48        double* melting=NULL;
     49       
    4850
    4951        /*diverse: */
    5052        int     numberofnodes;
    5153        int     analysis_type;
    52         int     sub_analysis_type;
    5354        int     count;
    5455
    5556        parameters->FindParam((void*)&analysis_type,"analysis_type");
    56         parameters->FindParam((void*)&sub_analysis_type,"sub_analysis_type");
    5757        count=parameters->Size();
    5858               
    … …  
    6363
    6464       
    65         if (   (analysis_type==ControlAnalysisEnum()) ||  (analysis_type==DiagnosticAnalysisEnum()) || (analysis_type==ThermalAnalysisEnum())){
     65        if (   (analysis_type==ControlAnalysisEnum()) || 
     66                   (analysis_type==DiagnosticAnalysisEnum()) ||
     67                   (analysis_type==ThermalAnalysisEnum()) ||
     68                   (analysis_type==PrognosticAnalysisEnum())
     69                   ){
    6670
    6771                parameters->FindParam((void*)&vx,"vx");
    … …  
    170174                param=(Param*)parameters->FindParamObject("pressure");
    171175                parameters->DeleteObject((Object*)param);
    172 
    173                 if(sub_analysis_type==TransientAnalysisEnum()){
    174 
    175                         parameters->FindParam((void*)&temperature,"temperature");
    176                         parameters->FindParam((void*)&melting,"melting");
    177 
    178                         /*Now, from temperature and melting, build t_g and m_g, correctly partitioned: */
    179                         t_g=(double*)xcalloc(numberofnodes,sizeof(double));
    180                         m_g=(double*)xcalloc(numberofnodes,sizeof(double));
    181 
    182                         for(i=0;i<numberofnodes;i++){
    183                                 t_g[(int)(partition[i])]=temperature[i];
    184                                 m_g[(int)(partition[i])]=melting[i];
    185                         }
    186 
    187                         /*Now, create new parameter: */
    188                         count++;
    189                         param= new Param(count,"t_g",DOUBLEVEC);
    190                         param->SetDoubleVec(t_g,numberofnodes);
    191                         parameters->AddObject(param);
    192 
    193                         count++;
    194                         param= new Param(count,"m_g",DOUBLEVEC);
    195                         param->SetDoubleVec(m_g,numberofnodes);
    196                         parameters->AddObject(param);
    197 
    198                         /*erase old parameter: */
    199                         param=(Param*)parameters->FindParamObject("temperature");
    200                         parameters->DeleteObject((Object*)param);
    201 
    202                         param=(Param*)parameters->FindParamObject("melting");
    203                         parameters->DeleteObject((Object*)param);
    204                 }
    205         }
    206 
     176        }
     177
     178        if(analysis_type==PrognosticAnalysisEnum()){
     179
     180                parameters->FindParam((void*)&melting,"melting");
     181               
     182                /*Now, from melting, build m_g, correctly partitioned: */
     183                m_g=(double*)xcalloc(numberofnodes,sizeof(double));
     184
     185                for(i=0;i<numberofnodes;i++){
     186                        m_g[(int)(partition[i])]= melting[i]; 
     187                }
     188
     189                /*Now, create new parameter: */
     190                count++;
     191                param= new Param(count,"m_g",DOUBLEVEC);
     192                param->SetDoubleVec(m_g,numberofnodes);
     193                parameters->AddObject(param);
     194
     195                /*erase old parameter: */
     196                param=(Param*)parameters->FindParamObject("melting");
     197                parameters->DeleteObject((Object*)param);
     198
     199               
     200               
     201                parameters->FindParam((void*)&thickness,"thickness");
     202               
     203                /*Now, from thickness, build h_g, correctly partitioned: */
     204                h_g=(double*)xcalloc(numberofnodes,sizeof(double));
     205
     206                for(i=0;i<numberofnodes;i++){
     207                        h_g[(int)(partition[i])]= thickness[i]; 
     208                }
     209
     210                /*Now, create new parameter: */
     211                count++;
     212                param= new Param(count,"h_g",DOUBLEVEC);
     213                param->SetDoubleVec(h_g,numberofnodes);
     214                parameters->AddObject(param);
     215
     216                /*erase old parameter: */
     217                param=(Param*)parameters->FindParamObject("thickness");
     218                parameters->DeleteObject((Object*)param);
     219
     220
     221
     222
     223                parameters->FindParam((void*)&accumulation,"accumulation");
     224               
     225                /*Now, from accumulation, build a_g, correctly partitioned: */
     226                a_g=(double*)xcalloc(numberofnodes,sizeof(double));
     227
     228                for(i=0;i<numberofnodes;i++){
     229                        a_g[(int)(partition[i])]= accumulation[i]; 
     230                }
     231
     232                /*Now, create new parameter: */
     233                count++;
     234                param= new Param(count,"a_g",DOUBLEVEC);
     235                param->SetDoubleVec(a_g,numberofnodes);
     236                parameters->AddObject(param);
     237
     238                /*erase old parameter: */
     239                param=(Param*)parameters->FindParamObject("accumulation");
     240                parameters->DeleteObject((Object*)param);
     241               
     242
     243
     244        }
    207245        xfree((void**)&partition);
    208246       
    … …  
    215253        xfree((void**)&vy_obs);
    216254        xfree((void**)&pressure);
    217         xfree((void**)&temperature);
    218         xfree((void**)&melting);
    219255        xfree((void**)&control_parameter);
    220256        xfree((void**)&u_g_obs);
    221257        xfree((void**)&p_g);
    222         xfree((void**)&t_g);
     258        xfree((void**)&a_g);
     259        xfree((void**)&h_g);
    223260        xfree((void**)&m_g);
    224261

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

    r358 r586  
    4545#include "./ControlConstrainx/ControlConstrainx.h"
    4646#include "./VelocityExtrudex/VelocityExtrudex.h"
     47#include "./VelocityDepthAveragex/VelocityDepthAveragex.h"
    4748#include "./GriddataMeshToGridx/GriddataMeshToGridx.h"
    4849#include "./ComputePressurex/ComputePressurex.h"
    4950#include "./SlopeExtrudex/SlopeExtrudex.h"
     51#include "./ThicknessExtrudex/ThicknessExtrudex.h"
    5052
    5153

  • issm/trunk/src/c/objects/Param.cpp

    r209 r586  
    1919               
    2020Param::Param(){
     21        stringarray=NULL;
    2122        doublevec=NULL;
    2223        doublemat=NULL;
    … …  
    3435        strcpy(name,param_name);
    3536        type=param_type;
    36         if ((param_type!=STRING) & (param_type!=INTEGER) & (param_type!=DOUBLE) &
     37        if ((param_type!=STRING) & (param_type!=STRINGARRAY) & (param_type!=INTEGER) & (param_type!=DOUBLE) &
    3738                (param_type!=DOUBLEVEC) & (param_type!=DOUBLEMAT) & (param_type!=PETSCVEC) & (param_type!=PETSCMAT)
    3839                ){
    3940                throw ErrorException(__FUNCT__,exprintf("%s%i"," unknow parameter type ",param_type));
    4041        }
     42        stringarray=NULL;
    4143        doublevec=NULL;
    4244        doublemat=NULL;
    … …  
    4749
    4850Param::~Param(){
     51       
     52        int i;
    4953
    5054        switch(type){
    … …  
    5256                case STRING:
    5357                        break;
     58
     59                case STRINGARRAY:
     60
     61                        for(i=0;i<M;i++){
     62                                char* descriptor=stringarray[i];
     63                                xfree((void**)&descriptor);
     64                        }
     65                        xfree((void**)&stringarray);
     66
     67                        break;
    5468       
    5569                case INTEGER:
    5670                        break;
     71
    5772       
    5873                case DOUBLE:
    … …  
    94109                        printf("   string value: %s\n",string);
    95110                        break;
     111                       
     112                case  STRINGARRAY:
     113                        printf("   string array: %i strings\n",M);
     114                        for(i=0;i<M;i++){
     115                                printf("      %i: %s\n",i,stringarray[i]);
     116                        }
    96117       
    97118                case INTEGER:
    … …  
    142163        double* serial_vec=NULL;
    143164        double* serial_mat=NULL;
     165        int i;
    144166
    145167        /*recover marshalled_dataset: */
    … …  
    161183                        memcpy(marshalled_dataset,&string,sizeof(string));marshalled_dataset+=sizeof(string);
    162184                        break;
     185
     186                case STRINGARRAY:
     187                        memcpy(marshalled_dataset,&M,sizeof(M));marshalled_dataset+=sizeof(M);
     188                        for(i=0;i<M;i++){
     189                                int size=(strlen(stringarray[i])+1)*sizeof(char);
     190                                memcpy(marshalled_dataset,&size,sizeof(size));marshalled_dataset+=sizeof(size);
     191                                memcpy(marshalled_dataset,&stringarray[i],size);marshalled_dataset+=size;
     192                        }
     193                        break;
     194
    163195                case INTEGER:
    164196                        memcpy(marshalled_dataset,&integer,sizeof(integer));marshalled_dataset+=sizeof(integer);
    … …  
    228260
    229261        int size;
     262        int i;
    230263
    231264        size=sizeof(id)+
    … …  
    238271                        size+=sizeof(string);
    239272                        break;
     273
     274                case STRINGARRAY:
     275                        size+=sizeof(M);
     276                        for(i=0;i<M;i++){
     277                                size+=sizeof(integer);
     278                                size+=(strlen(stringarray[i])+1)*sizeof(char);
     279                        }
     280                        break;
     281
    240282                case INTEGER:
    241283                        size+= sizeof(integer);
    … …  
    297339                case STRING:
    298340                        memcpy(&string,marshalled_dataset,sizeof(string));marshalled_dataset+=sizeof(string);
     341                        break;
     342               
     343                case STRINGARRAY:
     344                        memcpy(&M,marshalled_dataset,sizeof(M));marshalled_dataset+=sizeof(M);
     345                        if(M){
     346                                stringarray=(char**)xmalloc(M*sizeof(char*));
     347                                for(i=0;i<M;i++){
     348                                        int size;
     349                                        memcpy(&size,marshalled_dataset,sizeof(integer));marshalled_dataset+=sizeof(integer);
     350                                        memcpy(&stringarray[i],marshalled_dataset,size);marshalled_dataset+=size;
     351                                }
     352                        }
    299353                        break;
    300354
    … …  
    392446void  Param::GetParameterValue(void* pvalue){
    393447
     448        int i;
     449
    394450        if (type==STRING){
    395451                char** pstring=(char**)pvalue; //a little dangerous, but hey!
    … …  
    398454                strcpy(outstring,string);
    399455                *pstring=outstring;
     456        }
     457        else if (type==STRINGARRAY){
     458                char*** pstringarray=(char***)pvalue; //a little dangerous, but hey!
     459                char**  outstringarray=NULL;
     460                outstringarray=(char**)xmalloc(M*sizeof(char*));
     461                for(i=0;i<M;i++){
     462                        char* outstring=(char*)xmalloc((strlen(stringarray[i])+1)*sizeof(char));
     463                        strcpy(outstring,stringarray[i]);
     464                        outstringarray[i]=outstring;
     465                }
     466                *pstringarray=outstringarray;
    400467        }
    401468        else if(type==INTEGER){
    … …  
    488555}
    489556
     557#undef __FUNCT__
     558#define __FUNCT__ "SetStringArray"
     559void  Param::SetStringArray(char** value,int size){
     560       
     561        int i;
     562        if (type!=STRINGARRAY) throw ErrorException(__FUNCT__,exprintf("%s%i"," trying to set string for type",type));
     563        M=size;
     564        stringarray=(char**)xmalloc(M*sizeof(char*));
     565        for(i=0;i<M;i++){
     566                char* instring=(char*)xmalloc((strlen(value[i])+1)*sizeof(char));
     567                strcpy(instring,value[i]);
     568                stringarray[i]=instring;
     569        }
     570}
     571
     572
    490573int   Param::GetM(){
    491574        return M;

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

    r202 r586  
    1111#define PARAMSTRING 200 //max string length
    1212
    13 enum param_type { STRING, INTEGER, DOUBLE, DOUBLEVEC, DOUBLEMAT, PETSCVEC, PETSCMAT };
     13enum param_type { STRING, STRINGARRAY, INTEGER, DOUBLE, DOUBLEVEC, DOUBLEMAT, PETSCVEC, PETSCMAT };
    1414
    1515class Param: public Object{
    … …  
    2323                double ddouble;
    2424                char  string[PARAMSTRING];
     25                char** stringarray;
    2526                double* doublevec;
    2627                double* doublemat;
    … …  
    4950                void  SetInteger(int value);
    5051                void  SetString(char* value);
     52                void  SetStringArray(char** value,int size);
    5153                void  GetParameterValue(void* pvalue);
    5254               

  • issm/trunk/src/c/objects/Penta.cpp

    r545 r586  
    1 /*!\file Penta.c
     1/*!\file Penta.cpp
    22 * \brief: implementation of the Penta object
    33 */
    … …  
    2020}
    2121Penta::Penta( int penta_id, int penta_mid, int penta_mparid, int penta_node_ids[6], double penta_h[6], double penta_s[6], double penta_b[6], double penta_k[6], int penta_friction_type,
    22                         double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel,
    23                         int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6],
    24                         int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_stokesreconditioning){
    25 
     22                                double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface, double penta_meanvel,double penta_epsvel,
     23                                int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6],
     24                                int penta_artdiff, int penta_thermal_steadystate,double penta_viscosity_overshoot,double penta_stokesreconditioning){
     25       
    2626        int i;
    2727
    … …  
    8383        printf("   b=[%i,%i,%i,%i,%i,%i]\n",b[0],b[1],b[2],b[3],b[4],b[5]);
    8484        printf("   k=[%i,%i,%i,%i,%i,%i]\n",k[0],k[1],k[2],k[3],k[4],k[5]);
    85 
     85       
    8686        printf("   friction_type: %i\n",friction_type);
    8787        printf("   p: %g\n",p);
    … …  
    9393        printf("   epsvel: %g\n",epsvel);
    9494        printf("   collapse: %i\n",collapse);
    95 
     95       
    9696        printf("   melting=[%i,%i,%i,%i,%i,%i]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]);
    9797        printf("   accumulation=[%i,%i,%i,%i,%i,%i]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]);
    … …  
    114114        /*get enum type of Penta: */
    115115        enum_type=PentaEnum();
    116 
     116       
    117117        /*marshall enum: */
    118118        memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
    119 
     119       
    120120        /*marshall Penta data: */
    121121        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
    … …  
    149149        memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
    150150        memcpy(marshalled_dataset,&stokesreconditioning,sizeof(stokesreconditioning));marshalled_dataset+=sizeof(stokesreconditioning);
    151 
     151       
    152152        *pmarshalled_dataset=marshalled_dataset;
    153153        return;
    154154}
    155 
     155               
    156156int   Penta::MarshallSize(){
    157157
    158158        return sizeof(id)+
    159           sizeof(mid)+
    160           sizeof(mparid)+
    161           sizeof(node_ids)+
    162           sizeof(nodes)+
    163           sizeof(node_offsets)+
    164           sizeof(matice)+
    165           sizeof(matice_offset)+
    166           sizeof(matpar)+
    167           sizeof(matpar_offset)+
    168           sizeof(h)+
    169           sizeof(s)+
    170           sizeof(b)+
    171           sizeof(k)+
    172           sizeof(friction_type)+
    173           sizeof(p)+
    174           sizeof(q)+
    175           sizeof(shelf)+
    176           sizeof(onbed)+
    177           sizeof(onsurface)+
    178           sizeof(meanvel)+
    179           sizeof(epsvel)+
    180           sizeof(collapse)+
    181           sizeof(melting)+
    182           sizeof(accumulation)+
    183           sizeof(geothermalflux)+
    184           sizeof(artdiff)+
    185           sizeof(thermal_steadystate) +
    186           sizeof(viscosity_overshoot) +
    187           sizeof(stokesreconditioning) +
    188           sizeof(int); //sizeof(int) for enum type
     159                sizeof(mid)+
     160                sizeof(mparid)+
     161                sizeof(node_ids)+
     162                sizeof(nodes)+
     163                sizeof(node_offsets)+
     164                sizeof(matice)+
     165                sizeof(matice_offset)+
     166                sizeof(matpar)+
     167                sizeof(matpar_offset)+
     168                sizeof(h)+
     169                sizeof(s)+
     170                sizeof(b)+
     171                sizeof(k)+
     172                sizeof(friction_type)+
     173                sizeof(p)+
     174                sizeof(q)+
     175                sizeof(shelf)+
     176                sizeof(onbed)+
     177                sizeof(onsurface)+
     178                sizeof(meanvel)+
     179                sizeof(epsvel)+
     180                sizeof(collapse)+
     181                sizeof(melting)+
     182                sizeof(accumulation)+
     183                sizeof(geothermalflux)+
     184                sizeof(artdiff)+
     185                sizeof(thermal_steadystate) +
     186                sizeof(viscosity_overshoot) +
     187                sizeof(stokesreconditioning) +
     188                sizeof(int); //sizeof(int) for enum type
    189189}
    190190
    … …  
    262262
    263263        int i;
    264 
     264       
    265265        DataSet* loadsin=NULL;
    266266        DataSet* nodesin=NULL;
    … …  
    274274        /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
    275275        ResolvePointers((Object**)nodes,node_ids,node_offsets,6,nodesin);
    276 
     276       
    277277        /*Same for materials: */
    278278        ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
    … …  
    293293        }
    294294        else if (analysis_type==DiagnosticAnalysisEnum()){
    295 
     295       
    296296                if (sub_analysis_type==HorizAnalysisEnum()){
    297 
     297               
    298298                        CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
    299299                }
    300300                else if (sub_analysis_type==VertAnalysisEnum()){
    301 
     301               
    302302                        CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
    303303                }
    304304                else if (sub_analysis_type==StokesAnalysisEnum()){
    305 
     305               
    306306                        CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type);
    307 
     307               
    308308                }
    309309                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
    310310        }
    311311        else if (analysis_type==SlopeComputeAnalysisEnum()){
    312 
     312               
    313313                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
    314314        }
    315315        else if (analysis_type==ThermalAnalysisEnum()){
    316 
     316               
    317317                CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type);
    318318        }
    319319        else if (analysis_type==MeltingAnalysisEnum()){
    320 
     320               
    321321                CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type);
    322322        }
    … …  
    342342        int          doflist[numdof];
    343343        int          numberofdofspernode;
    344 
    345 
     344       
     345       
    346346        /* 3d gaussian points: */
    347347        int     num_gauss,ig;
    … …  
    359359        int     num_area_gauss;
    360360        double  gauss_weight;
    361 
     361       
    362362        /* 2d gaussian point: */
    363363        int     num_gauss2d;
    … …  
    391391        double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    392392        double Jdet;
    393 
     393       
    394394        /*slope: */
    395395        double  slope[2]={0.0,0.0};
    … …  
    454454                }
    455455
    456 #ifdef _DEBUGELEMENTS_
     456                #ifdef _DEBUGELEMENTS_
    457457                if(my_rank==RANK && id==ELID){
    458458                        printf("El id %i Rank %i PentaElement input list before gaussian loop: \n",ELID,RANK);
    … …  
    471471                        printf("   temperature_average [%g %g %g %g %g %g]\n",temperature_average_list[0],temperature_average_list[1],temperature_average_list[2],temperature_average_list[3],temperature_average_list[4],temperature_average_list[5]);
    472472                }
    473 #endif
     473                #endif
    474474
    475475                /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    … …  
    484484                GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    485485
    486 #ifdef _DEBUGGAUSS_
     486                #ifdef _DEBUGGAUSS_
    487487                if(my_rank==RANK && id==ELID){
    488488                        printf("El id %i Rank %i PentaElement gauss points\n",ELID,RANK);
    … …  
    494494                        }       
    495495                }
    496 #endif
     496                #endif
    497497                /* Start  looping on the number of gaussian points: */
    498498                for (ig1=0; ig1<num_area_gauss; ig1++){
    499499                        for (ig2=0; ig2<num_vert_gauss; ig2++){
    500 
     500                       
    501501                                /*Pick up the gaussian point: */
    502502                                gauss_weight1=*(area_gauss_weights+ig1);
    503503                                gauss_weight2=*(vert_gauss_weights+ig2);
    504504                                gauss_weight=gauss_weight1*gauss_weight2;
    505 
    506 
     505                               
     506                               
    507507                                gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    508508                                gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    … …  
    514514                                GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
    515515                                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
    516 
     516                       
    517517                                /*Get viscosity: */
    518518                                matice->GetViscosity3d(&viscosity, &epsilon[0]);
    … …  
    525525                                /* Get Jacobian determinant: */
    526526                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
    527 
     527       
    528528                                /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    529529                                  onto this scalar matrix, so that we win some computational time: */
    530 
     530                               
    531531                                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    532532                                D_scalar=newviscosity*gauss_weight*Jdet;
    … …  
    534534                                        D[i][i]=D_scalar;
    535535                                }
    536 
     536               
    537537                                /*  Do the triple product tB*D*Bprime: */
    538538                                TripleMultiply( &B[0][0],5,numdof,1,
    539                                                         &D[0][0],5,5,0,
    540                                                         &Bprime[0][0],5,numdof,0,
    541                                                         &Ke_gg_gaussian[0][0],0);
     539                                                &D[0][0],5,5,0,
     540                                                &Bprime[0][0],5,numdof,0,
     541                                                &Ke_gg_gaussian[0][0],0);
    542542
    543543                                /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
    … …  
    549549                        } //for (ig2=0; ig2<num_vert_gauss; ig2++)
    550550                } //for (ig1=0; ig1<num_area_gauss; ig1++)
    551 
     551               
    552552
    553553                /*Add Ke_gg to global matrix Kgg: */
    554554                MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    555 
     555       
    556556                //Deal with 2d friction at the bedrock interface
    557557                if((onbed && !shelf)){
    … …  
    567567
    568568        }
    569 
    570 cleanup_and_return:
     569               
     570        cleanup_and_return:
    571571        xfree((void**)&first_gauss_area_coord);
    572572        xfree((void**)&second_gauss_area_coord);
    … …  
    627627        /*recover pointers: */
    628628        inputs=(ParameterInputs*)vinputs;
    629 
     629       
    630630
    631631        /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness
    … …  
    660660
    661661        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    662 #ifdef _DEBUG_
     662        #ifdef _DEBUG_
    663663        for (i=0;i<num_area_gauss;i++){
    664664                _printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
    … …  
    667667                _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
    668668        }
    669 #endif
    670 
     669        #endif
     670       
    671671        /* Start  looping on the number of gaussian points: */
    672672        for (ig1=0; ig1<num_area_gauss; ig1++){
    673673                for (ig2=0; ig2<num_vert_gauss; ig2++){
    674 
     674               
    675675                        /*Pick up the gaussian point: */
    676676                        gauss_weight1=*(area_gauss_weights+ig1);
    677677                        gauss_weight2=*(vert_gauss_weights+ig2);
    678678                        gauss_weight=gauss_weight1*gauss_weight2;
    679 
     679                       
    680680                        gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    681681                        gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    … …  
    693693                        /*  Do the triple product tB*D*Bprime: */
    694694                        TripleMultiply( &B[0][0],1,numgrids,1,
    695                                                 &DL_scalar,1,1,0,
    696                                                 &Bprime[0][0],1,numgrids,0,
    697                                                 &Ke_gg_gaussian[0][0],0);
     695                                        &DL_scalar,1,1,0,
     696                                        &Bprime[0][0],1,numgrids,0,
     697                                        &Ke_gg_gaussian[0][0],0);
    698698
    699699                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    … …  
    708708        /*Add Ke_gg to global matrix Kgg: */
    709709        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    710 
    711 cleanup_and_return:
     710       
     711        cleanup_and_return:
    712712        xfree((void**)&first_gauss_area_coord);
    713713        xfree((void**)&second_gauss_area_coord);
    … …  
    740740        double         gravity,rho_ice,rho_water;
    741741
    742 
     742       
    743743        /*Collapsed formulation: */
    744744        Tria*  tria=NULL;
    … …  
    746746        /*Grid data: */
    747747        double        xyz_list[numgrids][3];
    748 
     748       
    749749        /*parameters: */
    750750        double             xyz_list_tria[numgrids2d][3];
    … …  
    770770        double     DLStokes[14][14]={0.0};
    771771        double     tLDStokes[numdof2d][14];
    772 
     772       
    773773        /* gaussian points: */
    774774        int     num_area_gauss;
    … …  
    792792        double  alpha2_list[numgrids2d];
    793793        double  alpha2_gauss;
    794 
     794       
    795795        ParameterInputs* inputs=NULL;
    796796
    … …  
    804804                }
    805805        }
    806 
     806       
    807807        /*recovre material parameters: */
    808808        rho_water=matpar->GetRhoWater();
    … …  
    818818
    819819        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    820                 get tria gaussian points as well as segment gaussian points. For tria gaussian
    821                 points, order of integration is 2, because we need to integrate the product tB*D*B'
    822                 which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    823                 points, same deal, which yields 3 gaussian points.*/
     820           get tria gaussian points as well as segment gaussian points. For tria gaussian
     821           points, order of integration is 2, because we need to integrate the product tB*D*B'
     822           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     823           points, same deal, which yields 3 gaussian points.*/
    824824
    825825        area_order=5;
    … …  
    845845                        /*Compute strain rate: */
    846846                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    847 
     847               
    848848                        /*Get viscosity: */
    849849                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    … …  
    883883                /*Build alpha2_list used by drag stiffness matrix*/
    884884                Friction* friction=NewFriction();
    885 
     885               
    886886                /*Initialize all fields: */
    887887                friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    888888                strcpy(friction->element_type,"2d");
    889 
     889               
    890890                friction->gravity=matpar->GetG();
    891891                friction->rho_ice=matpar->GetRhoIce();
    … …  
    910910                        }
    911911                }
    912 
     912               
    913913                GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
    914914
    … …  
    920920                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
    921921                        gauss_coord[3]=-1;
    922 
     922                       
    923923                        gauss_coord_tria[0]=*(first_gauss_area_coord+igarea);
    924924                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
    925925                        gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
    926 
     926       
    927927                        /*Get the Jacobian determinant */
    928928                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
    929 
     929                       
    930930                        /*Get L matrix if viscous basal drag present: */
    931931                        GetLStokes(&LStokes[0][0],  gauss_coord_tria);
    932932                        GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord);
    933 
     933                               
    934934                        /*Compute strain rate: */
    935935                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    936 
     936               
    937937                        /*Get viscosity at last iteration: */
    938938                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    939 
     939       
    940940                        /*Get normal vecyor to the bed */
    941941                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
    942 
     942                       
    943943                        bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
    944944                        bed_normal[1]=-surface_normal[1];
    … …  
    974974                }
    975975        } //if ( (onbed==1) && (shelf==0))
    976 
     976       
    977977        /*Reduce the matrix */
    978978        ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]);
    … …  
    986986        /*Add Ke_gg to global matrix Kgg: */
    987987        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
    988 
    989 cleanup_and_return:
     988       
     989        cleanup_and_return:
    990990
    991991        return;
    … …  
    998998        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    999999        if (analysis_type==ControlAnalysisEnum()){
    1000 
     1000               
    10011001                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
    10021002        }
    10031003        else if (analysis_type==DiagnosticAnalysisEnum()){
    1004 
     1004       
    10051005                if (sub_analysis_type==HorizAnalysisEnum()){
    10061006
    … …  
    10081008                }
    10091009                else if (sub_analysis_type==VertAnalysisEnum()){
    1010 
     1010               
    10111011                        CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
    10121012                }
    10131013                else if (sub_analysis_type==StokesAnalysisEnum()){
    1014 
     1014               
    10151015                        CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
    10161016                }
    … …  
    10181018        }
    10191019        else if (analysis_type==SlopeComputeAnalysisEnum()){
    1020 
     1020               
    10211021                CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
    10221022        }
    10231023        else if (analysis_type==ThermalAnalysisEnum()){
    1024 
     1024               
    10251025                CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type);
    10261026        }
    10271027        else if (analysis_type==MeltingAnalysisEnum()){
    1028 
     1028               
    10291029                CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type);
    10301030        }
    … …  
    10561056        inputs->Recover("melting",&melting[0],1,dofs,6,(void**)nodes);
    10571057        inputs->Recover("accumulation_param",&accumulation[0],1,dofs,6,(void**)nodes);
    1058 
     1058       
    10591059        //Update material if necessary
    10601060        if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){
    … …  
    10631063                matice->SetB(B_average);
    10641064        }
    1065 
     1065       
    10661066        if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){
    10671067                B_average=(B_list[0]+B_list[1]+B_list[2]+B_list[3]+B_list[4]+B_list[5])/6.0;
    … …  
    10881088        }
    10891089}
    1090 
     1090               
    10911091int Penta::GetOnBed(){
    10921092        return onbed;
    … …  
    10991099}
    11001100void          Penta::GetBedList(double* bed_list){
    1101 
     1101       
    11021102        int i;
    11031103        for(i=0;i<6;i++)bed_list[i]=b[i];
    … …  
    11151115        int i;
    11161116        Tria* tria=NULL;
    1117 
     1117       
    11181118        /*Bail out if this element if:
    11191119         * -> Non collapsed and not on the surface
    … …  
    11321132        }
    11331133        else{
    1134 
     1134               
    11351135                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    11361136                tria->Du(du_g,u_g,u_g_obs,inputs,analysis_type,sub_analysis_type);
    … …  
    11561156#define __FUNCT__ "Penta::GradjDrag"
    11571157void  Penta::GradjDrag(Vec grad_g,double* u_g,double* lambda_g,void* inputs,int analysis_type,int sub_analysis_type){
    1158 
    1159 
     1158       
     1159       
    11601160        Tria* tria=NULL;
    1161 
     1161       
    11621162        /*Bail out if this element does not touch the bedrock: */
    11631163        if (!onbed){
    … …  
    11651165        }
    11661166        else{
    1167 
     1167               
    11681168                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    11691169                tria->GradjDrag( grad_g,u_g,lambda_g,inputs,analysis_type,sub_analysis_type);
    … …  
    11781178        throw ErrorException(__FUNCT__," not supported yet!");
    11791179}
    1180 
     1180       
    11811181#undef __FUNCT__
    11821182#define __FUNCT__ "Penta::Misfit"
    11831183double Penta::Misfit(double* velocity,double* obs_velocity,void* inputs,int analysis_type,int sub_analysis_type){
    1184 
     1184       
    11851185        double J;
    11861186        Tria* tria=NULL;
    1187 
     1187       
    11881188        /*Bail out if this element if:
    11891189         * -> Non collapsed and not on the surface
    … …  
    12091209        }
    12101210}
    1211 
     1211               
    12121212#undef __FUNCT__
    12131213#define __FUNCT__ "Penta::SpawnTria"
    … …  
    12241224        double   tria_accumulation[3];
    12251225        double   tria_geothermalflux[3];
    1226 
     1226       
    12271227        /*configuration: */
    12281228        int tria_node_ids[3];
    … …  
    12491249        tria_melting[1]=melting[g1];
    12501250        tria_melting[2]=melting[g2];
    1251 
     1251       
    12521252        tria_accumulation[0]=accumulation[g0];
    12531253        tria_accumulation[1]=accumulation[g1];
    … …  
    12701270        tria_node_offsets[2]=node_offsets[g2];
    12711271
    1272         tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot);
     1272        tria= new Tria(id,mid,mparid,tria_node_ids,tria_h,tria_s,tria_b,tria_k, tria_melting, tria_accumulation, tria_geothermalflux,friction_type,p,q,shelf,meanvel,epsvel,viscosity_overshoot,artdiff);
    12731273
    12741274        tria->NodeConfiguration(tria_node_ids,tria_nodes,tria_node_offsets);
    … …  
    12861286        int doflist_per_node[MAXDOFSPERNODE];
    12871287        int numberofdofspernode;
    1288 
     1288       
    12891289        for(i=0;i<6;i++){
    12901290                nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
    … …  
    13201320        GetB(&B[0][0], xyz_list, gauss_l1l2l3l4);
    13211321
    1322 #ifdef _DEBUG_
     1322        #ifdef _DEBUG_
    13231323        _printf_("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
    13241324        _printf_("              [ %lf   %lf  \n",B[1][0],B[1][1]);
    … …  
    13261326        _printf_("              [ %lf   %lf  ]\n",B[3][0],B[3][1]);
    13271327        _printf_("              [ %lf   %lf  ]\n",B[4][0],B[4][1]);
    1328 
     1328       
    13291329        _printf_("B for grid2 : [ %lf   %lf  \n",B[0][2],B[0][3]);
    13301330        _printf_("              [ %lf   %lf  \n",B[1][2],B[1][3]);
    … …  
    13321332        _printf_("              [ %lf   %lf  ]\n",B[3][2],B[3][3]);
    13331333        _printf_("              [ %lf   %lf  ]\n",B[4][2],B[4][3]);
    1334 
     1334       
    13351335        _printf_("B for grid3 : [ %lf   %lf  \n", B[0][4],B[0][5]);
    13361336        _printf_("              [ %lf   %lf  \n", B[1][4],B[1][5]);
    … …  
    13381338        _printf_("              [ %lf   %lf  ]\n",B[3][4],B[3][5]);
    13391339        _printf_("              [ %lf   %lf  ]\n",B[4][4],B[4][5]);
    1340 
     1340       
    13411341        _printf_("B for grid4 : [ %lf   %lf  \n", B[0][6],B[0][7]);
    13421342        _printf_("              [ %lf   %lf  \n", B[1][6],B[1][7]);
    … …  
    13441344        _printf_("              [ %lf   %lf  ]\n",B[3][6],B[3][7]);
    13451345        _printf_("              [ %lf   %lf  ]\n",B[4][6],B[4][7]);
    1346 
     1346                               
    13471347        _printf_("B for grid5 : [ %lf   %lf  \n", B[0][8],B[0][9]);
    13481348        _printf_("              [ %lf   %lf  \n", B[1][8],B[1][9]);
    … …  
    13601360                _printf_("Velocity for grid %i %lf %lf\n",i,*(vxvy_list+2*i+0),*(vxvy_list+2*i+1));
    13611361        }
    1362 #endif
     1362        #endif
    13631363
    13641364        /*Multiply B by velocity, to get strain rate: */
    13651365        MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
    1366                                 velocity,NDOF2*numgrids,1,0,
    1367                                 epsilon,0);
     1366                                      velocity,NDOF2*numgrids,1,0,
     1367                                                  epsilon,0);
    13681368
    13691369}
    … …  
    13851385         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    13861386         */
    1387 
     1387       
    13881388        int i;
    13891389        const int numgrids=6;
    … …  
    13961396        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    13971397
    1398 #ifdef _DEBUG_
     1398        #ifdef _DEBUG_
    13991399        for (i=0;i<numgrids;i++){
    14001400                _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
    14011401        }
    1402 #endif
     1402        #endif
    14031403
    14041404        /*Build B: */
    … …  
    14151415                *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
    14161416                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    1417 
     1417               
    14181418                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    14191419                *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
    … …  
    14391439         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    14401440         */
    1441 
     1441       
    14421442        int i;
    14431443        const int NDOF3=3;
    … …  
    14501450        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    14511451
    1452 #ifdef _DEBUG_
     1452        #ifdef _DEBUG_
    14531453        for (i=0;i<numgrids;i++){
    14541454                _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
    14551455        }
    1456 #endif
     1456        #endif
    14571457
    14581458        /*Build BPrime: */
    … …  
    14691469                *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[2][i];
    14701470                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    1471 
     1471               
    14721472                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    14731473                *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[2][i];
    … …  
    14821482         * the determinant of it: */
    14831483        const int NDOF3=3;
    1484 
     1484       
    14851485        double J[NDOF3][NDOF3];
    14861486
    … …  
    14961496#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasic"
    14971497void Penta::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4){
    1498 
     1498       
    14991499        /*This routine returns the values of the nodal functions derivatives  (with respect to the basic coordinate system: */
    15001500
    1501 
     1501       
    15021502        int i;
    15031503        const int NDOF3=3;
    … …  
    15341534        const int NDOF3=3;
    15351535        int i,j;
    1536 
     1536       
    15371537        /*The Jacobian is constant over the element, discard the gaussian points.
    15381538         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    … …  
    15441544        double y1,y2,y3,y4,y5,y6;
    15451545        double z1,z2,z3,z4,z5,z6;
    1546 
     1546       
    15471547        double sqrt3=sqrt(3.0);
    1548 
     1548       
    15491549        /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
    15501550        A1=gauss_l1l2l3l4[0];
    … …  
    15621562        x5=*(xyz_list+3*4+0);
    15631563        x6=*(xyz_list+3*5+0);
    1564 
     1564       
    15651565        y1=*(xyz_list+3*0+1);
    15661566        y2=*(xyz_list+3*1+1);
    … …  
    15901590        *(J+NDOF3*2+2)=sqrt3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+1.0/4.0*(z1-z2-z4+z5)*xi+1.0/4.0*(-z1+z5-z2+z4);
    15911591
    1592 #ifdef _DEBUG_
     1592        #ifdef _DEBUG_
    15931593        for(i=0;i<3;i++){
    15941594                for (j=0;j<3;j++){
    … …  
    15971597                printf("\n");
    15981598        }
    1599 #endif
     1599        #endif
    16001600}
    16011601
    … …  
    16031603#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParams"
    16041604void Penta::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4){
    1605 
     1605       
    16061606        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    16071607         * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
    … …  
    16101610        double A1,A2,A3,z;
    16111611        double sqrt3=sqrt(3.0);
    1612 
     1612       
    16131613        A1=gauss_l1l2l3l4[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
    16141614        A2=gauss_l1l2l3l4[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
    … …  
    16771677        int          doflist[numdof];
    16781678        int          numberofdofspernode;
    1679 
     1679       
    16801680        /* parameters: */
    16811681        double  slope[NDOF2];
    … …  
    17521752
    17531753                GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    1754 #ifdef _DEBUG_
     1754                #ifdef _DEBUG_
    17551755                for (i=0;i<num_area_gauss;i++){
    17561756                        _printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
    … …  
    17591759                        _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
    17601760                }
    1761 #endif
    1762 
     1761                #endif
     1762       
    17631763                /* Start  looping on the number of gaussian points: */
    17641764                for (ig1=0; ig1<num_area_gauss; ig1++){
    17651765                        for (ig2=0; ig2<num_vert_gauss; ig2++){
    1766 
     1766                       
    17671767                                /*Pick up the gaussian point: */
    17681768                                gauss_weight1=*(area_gauss_weights+ig1);
    17691769                                gauss_weight2=*(vert_gauss_weights+ig2);
    17701770                                gauss_weight=gauss_weight1*gauss_weight2;
    1771 
     1771                               
    17721772                                gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    17731773                                gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    17741774                                gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    17751775                                gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
    1776 
     1776               
    17771777                                /*Compute thickness at gaussian point: */
    17781778                                GetParameterValue(&thickness, &h[0],gauss_l1l2l3l4);
    … …  
    17831783                                /* Get Jacobian determinant: */
    17841784                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
    1785 
    1786                                 /*Get nodal functions: */
     1785               
     1786                                 /*Get nodal functions: */
    17871787                                GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
    17881788
    … …  
    18081808        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    18091809
    1810 cleanup_and_return:
     1810        cleanup_and_return:
    18111811        xfree((void**)&first_gauss_area_coord);
    18121812        xfree((void**)&second_gauss_area_coord);
    … …  
    18221822#define __FUNCT__ "Penta::CreateKMatrix"
    18231823void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
    1824 
     1824       
    18251825        const int numgrids=6;
    18261826        double l1l2l3l4l5l6[numgrids];
    … …  
    18341834#define __FUNCT__ "Penta::GetParameterDerivativeValue"
    18351835void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){
    1836 
     1836                               
    18371837        /*From grid values of parameter p (p_list[0], p_list[1], p_list[2], p_list[3], p_list[4] and p_list[4]), return parameter derivative value at gaussian point specified by gauss_l1l2l3l4:
    18381838         *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
    … …  
    18491849        /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
    18501850        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    1851 
     1851       
    18521852        *(p+0)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[0][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[0][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[0][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[0][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[0][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[0][5];
    1853         ;
     1853;
    18541854        *(p+1)=p_list[0]*dh1dh2dh3dh4dh5dh6_basic[1][0]+p_list[1]*dh1dh2dh3dh4dh5dh6_basic[1][1]+p_list[2]*dh1dh2dh3dh4dh5dh6_basic[1][2]+p_list[3]*dh1dh2dh3dh4dh5dh6_basic[1][3]+p_list[4]*dh1dh2dh3dh4dh5dh6_basic[1][4]+p_list[5]*dh1dh2dh3dh4dh5dh6_basic[1][5];
    18551855
    … …  
    18611861#define __FUNCT__ "Penta::GetNodalFunctions"
    18621862void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){
    1863 
     1863       
    18641864        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    18651865
    … …  
    18671867
    18681868        l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0;
    1869 
     1869       
    18701870        l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0;
    18711871
    … …  
    18731873
    18741874        l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0;
    1875 
     1875       
    18761876        l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0;
    18771877
    … …  
    18881888        int          nodedofs[2];
    18891889        int          numberofdofspernode;
    1890 
     1890       
    18911891        Node* node=NULL;
    18921892        int   i;
    18931893        double velocity[2];
    1894 
     1894               
    18951895
    18961896        if((collapse==1) && (onbed==1)){
    1897 
     1897                       
    18981898                GetDofList(&doflist[0],&numberofdofspernode);
    18991899
    … …  
    19021902                 * inserting the same velocity value into ug, until we reach the surface: */
    19031903                for(i=0;i<3;i++){
    1904 
     1904       
    19051905                        node=nodes[i]; //base nodes
    1906 
     1906       
    19071907                        /*get velocity for this base node: */
    19081908                        velocity[0]=ug_serial[doflist[numberofdofspernode*i+0]];
    … …  
    19261926
    19271927}
     1928
     1929#undef __FUNCT__
     1930#define __FUNCT__ "Penta::VelocityExtrudeAllElements"
     1931void  Penta::VelocityExtrudeAllElements(Vec ug,double* ug_serial){
     1932
     1933        /* node data: */
     1934        const int    numgrids=6;
     1935        const int    numdof=2*numgrids;
     1936        int          doflist[numdof];
     1937        int          nodedofs[2];
     1938        int          numberofdofspernode;
     1939       
     1940        Node* node=NULL;
     1941        int   i;
     1942        double velocity[2];
     1943               
     1944
     1945        if(onbed==1){
     1946                       
     1947                GetDofList(&doflist[0],&numberofdofspernode);
     1948
     1949                /*this penta is on the bed. For each node on the base of this penta,
     1950                 * we grab the velocity. Once we know the velocity, we follow the upper nodes,
     1951                 * inserting the same velocity value into ug, until we reach the surface: */
     1952                for(i=0;i<3;i++){
     1953       
     1954                        node=nodes[i]; //base nodes
     1955       
     1956                        /*get velocity for this base node: */
     1957                        velocity[0]=ug_serial[doflist[numberofdofspernode*i+0]];
     1958                        velocity[1]=ug_serial[doflist[numberofdofspernode*i+1]];
     1959
     1960                        //go througn all nodes which sit on top of this node, until we reach the surface,
     1961                        //and plug  velocity in ug
     1962                        for(;;){
     1963
     1964                                node->GetDofList(&nodedofs[0],&numberofdofspernode);
     1965                                VecSetValues(ug,1,&nodedofs[0],&velocity[0],INSERT_VALUES);
     1966                                VecSetValues(ug,1,&nodedofs[1],&velocity[1],INSERT_VALUES);
     1967
     1968                                if (node->IsOnSurface())break;
     1969                                /*get next node: */
     1970                                node=node->GetUpperNode();
     1971                        }
     1972                }
     1973
     1974        }
     1975
     1976}
     1977
     1978#undef __FUNCT__
     1979#define __FUNCT__ "Penta::ThicknessExtrude"
     1980void  Penta::ThicknessExtrude(Vec tg,double* tg_serial){
     1981
     1982        /* node data: */
     1983        const int    numgrids=6;
     1984        const int    numdof=1*numgrids;
     1985        int          doflist[numdof];
     1986        int          nodedofs;
     1987        int          numberofdofspernode;
     1988       
     1989        Node* node=NULL;
     1990        int   i;
     1991        double thickness;
     1992               
     1993
     1994        if(onbed==1){
     1995                       
     1996                GetDofList(&doflist[0],&numberofdofspernode);
     1997
     1998                /*this penta is on the bed. For each node on the base of this penta,
     1999                 * we grab the thickness. Once we know the thickness, we follow the upper nodes,
     2000                 * inserting the same thickness value into tg, until we reach the surface: */
     2001                for(i=0;i<3;i++){
     2002       
     2003                        node=nodes[i]; //base nodes
     2004       
     2005                        /*get velocity for this base node: */
     2006                        thickness=tg_serial[doflist[numberofdofspernode*i+0]];
     2007
     2008                        //go through all nodes which sit on top of this node, until we reach the surface,
     2009                        //and pltg  velocity in tg
     2010                        for(;;){
     2011
     2012                                node->GetDofList(&nodedofs,&numberofdofspernode);
     2013                                VecSetValues(tg,1,&nodedofs,&thickness,INSERT_VALUES);
     2014
     2015                                if (node->IsOnSurface())break;
     2016                                /*get next node: */
     2017                                node=node->GetUpperNode();
     2018                        }
     2019                }
     2020
     2021        }
     2022
     2023}
     2024
     2025#undef __FUNCT__
     2026#define __FUNCT__ "Penta::VelocityDepthAverageAtBase"
     2027void  Penta::VelocityDepthAverageAtBase(Vec ug,double* ug_serial){
     2028
     2029        /* node data: */
     2030        const int    numgrids=6;
     2031        const int    numdof=2*numgrids;
     2032        int          doflist[numdof];
     2033        int          nodedofs[2];
     2034        int          numberofdofspernode;
     2035        int          dofx,dofy;
     2036       
     2037        Node* node=NULL;
     2038        Node* upper_node=NULL;
     2039        int   i;
     2040        double base_velocity[2];
     2041        double velocity2[2];
     2042        double velocity1[2];
     2043        double velocity_average[2];
     2044        double sum[2];
     2045        double z1,z2,dz;
     2046        double thickness;
     2047               
     2048        if(onbed==1){
     2049                       
     2050                GetDofList(&doflist[0],&numberofdofspernode);
     2051
     2052                /*this penta is on the bed. For each node on the base of this penta, we are going to
     2053                 * follow the upper nodes until we reach the surface. At each upper node, we'll grab the
     2054                 * velocity for this node, and add it to ug. We'll finish by
     2055                 * we grab the velocity. Once we know the velocity, we follow the upper nodes,
     2056                 * inserting the same velocity value into ug, until we reach the surface: */
     2057                for(i=0;i<3;i++){
     2058
     2059                        //get thickness for this grid
     2060                        thickness=h[i]; //pick up from the penta element itself.
     2061       
     2062                        node=nodes[i]; //base nodes
     2063       
     2064                        /*get dofs for this base node velocity: */
     2065                        dofx=doflist[numberofdofspernode*i+0];
     2066                        dofy=doflist[numberofdofspernode*i+1];
     2067
     2068                        /*first thing, cancel velocity on the base nodes, so that we start with the value 0: */
     2069                        base_velocity[0]=-ug_serial[doflist[numberofdofspernode*i+0]]; //- to cancel
     2070                        base_velocity[1]=-ug_serial[doflist[numberofdofspernode*i+1]];
     2071               
     2072                        VecSetValues(ug,1,&dofx,&base_velocity[0],ADD_VALUES);
     2073                        VecSetValues(ug,1,&dofy,&base_velocity[1],ADD_VALUES);
     2074
     2075                        //go through all nodes which sit on top of this node, until we reach the surface,
     2076                        //and plug  deltaV*deltaH/H at base of ug: */
     2077                       
     2078                        for(;;){
     2079                               
     2080                                if (node->IsOnSurface())break;
     2081
     2082                                node->GetDofList(&nodedofs[0],&numberofdofspernode);
     2083                               
     2084                                velocity1[0]=ug_serial[nodedofs[numberofdofspernode*i+0]];
     2085                                velocity1[1]=ug_serial[nodedofs[numberofdofspernode*i+1]];
     2086                                z1=node->GetZ();
     2087
     2088                                upper_node=node->GetUpperNode();
     2089                                upper_node->GetDofList(&nodedofs[0],&numberofdofspernode);
     2090                       
     2091                                velocity2[0]=ug_serial[nodedofs[numberofdofspernode*i+0]];
     2092                                velocity2[1]=ug_serial[nodedofs[numberofdofspernode*i+1]];
     2093                                z2=upper_node->GetZ();
     2094
     2095                                dz=(z2-z1);
     2096                                velocity_average[0]=(velocity1[0]+velocity2[0])/2.0;
     2097                                velocity_average[1]=(velocity1[1]+velocity2[1])/2.0;
     2098
     2099                                sum[0]=velocity_average[0]*dz/thickness;
     2100                                sum[1]=velocity_average[1]*dz/thickness;
     2101
     2102                                /*Plug velocity_average*deltaH/H into base of ug: */
     2103                                VecSetValues(ug,1,&dofx,&sum[0],ADD_VALUES);
     2104                                VecSetValues(ug,1,&dofy,&sum[1],ADD_VALUES);
     2105
     2106                                /*get next node: */
     2107                                node=node->GetUpperNode();
     2108                        }
     2109                }
     2110
     2111        }
     2112
     2113}
     2114
    19282115
    19292116#undef __FUNCT__
    … …  
    19382125        int          nodedof;
    19392126        int          numberofdofspernode;
    1940 
     2127       
    19412128        Node* node=NULL;
    19422129        int   i;
    19432130        double slope;
    1944 
     2131               
    19452132
    19462133        if(onbed==1){
    1947 
     2134                       
    19482135                GetDofList(&doflist[0],&numberofdofspernode);
    19492136
    … …  
    19522139                /*For each node on the base of this penta,  we grab the slope. Once we know the slope, we follow the upper nodes,
    19532140                 * inserting the same slope value into sg, until we reach the surface: */
    1954 
     2141               
    19552142                for(i=0;i<3;i++){
    1956 
     2143       
    19572144                        node=nodes[i]; //base nodes
    1958 
     2145       
    19592146                        /*get velocity for this base node: */
    19602147                        slope=sg_serial[doflist[i]];
    … …  
    19842171
    19852172        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
    1986                 where hi is the interpolation function for grid i.*/
     2173        where hi is the interpolation function for grid i.*/
    19872174
    19882175        int i;
    … …  
    19942181        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    19952182
    1996 #ifdef _DEBUG_
     2183        #ifdef _DEBUG_
    19972184        for (i=0;i<numgrids;i++){
    19982185                _printf_("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf\n",i,dh1dh2dh3dh4dh5dh6_basic[0][i],dh1dh2dh3dh4dh5dh6_basic[1][i],dh1dh2dh3dh4dh5dh6_basic[2][i]);
    19992186        }
    2000 #endif
     2187        #endif
    20012188
    20022189        /*Build B: */
    … …  
    20042191                B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
    20052192        }
    2006 
     2193       
    20072194}
    20082195
    … …  
    20152202
    20162203        int i;
    2017 
     2204                               
    20182205        GetNodalFunctions(B, gauss_l1l2l3l4);
    20192206
    … …  
    20342221        int          doflist[numdof];
    20352222        int          numberofdofspernode;
    2036 
     2223       
    20372224        /* gaussian points: */
    20382225        int     num_gauss,ig;
    … …  
    20942281        /* recover input parameters: */
    20952282        if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
    2096         inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
     2283            inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
    20972284
    20982285        /*Get gaussian points and weights :*/
    … …  
    21012288
    21022289        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, &fourth_gauss_vert_coord,&vert_gauss_weights,order_area_gauss,num_vert_gauss);
    2103 #ifdef _DEBUG_
     2290        #ifdef _DEBUG_
    21042291        for (i=0;i<num_area_gauss;i++){
    21052292                _printf_("Area Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(area_gauss_weights+i));
    … …  
    21082295                _printf_("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
    21092296        }
    2110 #endif
     2297        #endif
    21112298
    21122299        /* Start  looping on the number of gaussian points: */
    21132300        for (ig1=0; ig1<num_area_gauss; ig1++){
    21142301                for (ig2=0; ig2<num_vert_gauss; ig2++){
    2115 
     2302               
    21162303                        /*Pick up the gaussian point: */
    21172304                        gauss_weight1=*(area_gauss_weights+ig1);
    21182305                        gauss_weight2=*(vert_gauss_weights+ig2);
    21192306                        gauss_weight=gauss_weight1*gauss_weight2;
    2120 
     2307                       
    21212308                        gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    21222309                        gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    21232310                        gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    21242311                        gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
    2125 
     2312       
    21262313                        /*Get velocity derivative, with respect to x and y: */
    21272314                        GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
    … …  
    21292316                        dudx=du[0];
    21302317                        dvdy=dv[1];
    2131 
     2318                       
    21322319
    21332320                        /* Get Jacobian determinant: */
    21342321                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
    2135 #ifdef _DEBUG_
     2322                        #ifdef _DEBUG_
    21362323                        _printf_("Element id %i Jacobian determinant: %lf\n",PentaElementGetID(this),Jdet);
    2137 #endif
    2138 
    2139                         /*Get nodal functions: */
     2324                        #endif
     2325               
     2326                         /*Get nodal functions: */
    21402327                        GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
    21412328
    … …  
    21542341        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    21552342
    2156 cleanup_and_return:
     2343        cleanup_and_return:
    21572344        xfree((void**)&first_gauss_area_coord);
    21582345        xfree((void**)&second_gauss_area_coord);
    … …  
    21752362        double rho_ice,g;
    21762363        double       xyz_list[numgrids][3];
    2177 
     2364               
    21782365        /*Get node data: */
    21792366        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    2180 
     2367       
    21812368        /*pressure is lithostatic: */
    21822369        //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab
    … …  
    21912378                pressure[i]=rho_ice*g*(s[i]-xyz_list[i][2]);
    21922379        }
    2193 
     2380       
    21942381        /*plug local pressure values into global pressure vector: */
    21952382        VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
    … …  
    22052392        /*Collapsed formulation: */
    22062393        Tria*  tria=NULL;
    2207 
     2394       
    22082395        /*Is this element on the bed? :*/
    22092396        if(!onbed)return;
    … …  
    22212408
    22222409void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    2223 
     2410       
    22242411        /*Collapsed formulation: */
    22252412        Tria*  tria=NULL;
    2226 
     2413       
    22272414        /*Is this element on the bed? :*/
    22282415        if(!onbed)return;
    … …  
    22382425#define __FUNCT__ "ReduceMatrixStokes"
    22392426void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
    2240 
     2427               
    22412428        int i,j;
    22422429
    … …  
    22922479        double det;
    22932480        int calculationdof=3;
    2294 
     2481       
    22952482        /*Take the matrix components: */
    22962483        a=*(Ke+calculationdof*0+0);
    … …  
    23052492
    23062493        det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g);
    2307 
     2494       
    23082495        *(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det;
    23092496        *(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det;
    … …  
    23922579         * For grid i, Bi can be expressed in the basic coordinate system
    23932580         * by:                          Bi_basic=[ dh/dx          0             0       0  ]
    2394          *                                      [   0           dh/dy           0       0  ]
    2395          *                                      [   0             0           dh/dy     0  ]
    2396          *                                      [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
    2397          *                                      [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
    2398          *                                      [   0          1/2*dh/dz    1/2*dh/dy   0  ]
    2399          *                                      [   0             0             0       h  ]
    2400          *                                      [ dh/dx         dh/dy         dh/dz     0  ]
    2401          *      where h is the interpolation function for grid i.
    2402          *      Same thing for Bb except the last column that does not exist.
    2403          */
    2404 
     2581                *                                       [   0           dh/dy           0       0  ]
     2582                *                                       [   0             0           dh/dy     0  ]
     2583                *                                       [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
     2584                *                                       [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
     2585                *                                       [   0          1/2*dh/dz    1/2*dh/dy   0  ]
     2586                *                                       [   0             0             0       h  ]
     2587                *                                       [ dh/dx         dh/dy         dh/dz     0  ]
     2588                *       where h is the interpolation function for grid i.
     2589                *       Same thing for Bb except the last column that does not exist.
     2590                */
     2591       
    24052592        int i;
    24062593        int calculationdof=3;
    … …  
    24162603
    24172604        GetNodalFunctions(l1l6, gauss_coord);
    2418 
    2419 #ifdef _DEBUG_
     2605       
     2606        #ifdef _DEBUG_
    24202607        for (i=0;i<7;i++){
    24212608                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i],dh1dh7_basic[2][i]);
    24222609        }
    24232610
    2424 #endif
     2611        #endif
    24252612
    24262613        /*Build B: */
    … …  
    24702657
    24712658        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    2472          *      For grid i, Bi' can be expressed in the basic coordinate system
    2473          *      by:
    2474          *                              Bi_basic'=[ dh/dx   0          0       0]
    2475          *                                       [   0      dh/dy      0       0]
    2476          *                                       [   0      0         dh/dz    0]
    2477          *                                       [  dh/dy   dh/dx      0       0]
    2478          *                                       [  dh/dz   0        dh/dx     0]
    2479          *                                       [   0      dh/dz    dh/dy     0]
    2480          *                                       [  dh/dx   dh/dy    dh/dz     0]
    2481          *                                       [   0      0          0       h]
    2482          *      where h is the interpolation function for grid i.
    2483          *
    2484          *      Same thing for the bubble fonction except that there is no fourth column
    2485          */
     2659        *       For grid i, Bi' can be expressed in the basic coordinate system
     2660        *       by:
     2661        *                               Bi_basic'=[ dh/dx   0          0       0]
     2662        *                                        [   0      dh/dy      0       0]
     2663        *                                        [   0      0         dh/dz    0]
     2664        *                                        [  dh/dy   dh/dx      0       0]
     2665        *                                        [  dh/dz   0        dh/dx     0]
     2666        *                                        [   0      dh/dz    dh/dy     0]
     2667        *                                        [  dh/dx   dh/dy    dh/dz     0]
     2668        *                                        [   0      0          0       h]
     2669        *       where h is the interpolation function for grid i.
     2670        *
     2671        *       Same thing for the bubble fonction except that there is no fourth column
     2672        */
    24862673
    24872674        int i;
    … …  
    24982685
    24992686        GetNodalFunctions(l1l6, gauss_coord);
    2500 
    2501 #ifdef _DEBUG_
     2687       
     2688        #ifdef _DEBUG_
    25022689        for (i=0;i<6;i++){
    25032690                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
    25042691        }
    25052692
    2506 #endif
     2693        #endif
    25072694
    25082695        /*B_primeuild B_prime: */
    … …  
    25512738
    25522739        /*
    2553          * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    2554          * For grid i, Li can be expressed in the basic coordinate system
    2555          * by:
    2556          *       Li_basic=[ h    0    0   0]
    2557          *               [ 0    h    0   0]
    2558          *               [ 0    0    h   0]
    2559          *               [ 0    0    h   0]
    2560          *               [ h    0    0   0]
    2561          *               [ 0    h    0   0]
    2562          *               [ h    0    0   0]
    2563          *               [ 0    h    0   0]
    2564          *               [ 0    0    h   0]
    2565          *               [ 0    0    h   0]
    2566          *               [ 0    0    h   0]
    2567          *               [ h    0    0   0]
    2568          *               [ 0    h    0   0]
    2569          *               [ 0    0    h   0]
    2570          * where h is the interpolation function for grid i.
    2571          */
     2740        * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     2741        * For grid i, Li can be expressed in the basic coordinate system
     2742        * by:
     2743        *       Li_basic=[ h    0    0   0]
     2744        *                [ 0    h    0   0]
     2745        *                [ 0    0    h   0]
     2746        *                [ 0    0    h   0]
     2747        *                [ h    0    0   0]
     2748        *                [ 0    h    0   0]
     2749        *                [ h    0    0   0]
     2750        *                [ 0    h    0   0]
     2751        *                [ 0    0    h   0]
     2752        *                [ 0    0    h   0]
     2753        *                [ 0    0    h   0]
     2754        *                [ h    0    0   0]
     2755        *                [ 0    h    0   0]
     2756        *                [ 0    0    h   0]
     2757        * where h is the interpolation function for grid i.
     2758        */
    25722759
    25732760        int i;
    … …  
    25822769        l1l2l3[1]=gauss_coord_tria[1];
    25832770        l1l2l3[2]=gauss_coord_tria[2];
    2584 
    2585 
    2586 #ifdef _DELUG_
     2771       
     2772
     2773        #ifdef _DELUG_
    25872774        for (i=0;i<3;i++){
    25882775                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
    25892776        }
    2590 #endif
     2777        #endif
    25912778
    25922779        /*Build LStokes: */
    … …  
    26482835                *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
    26492836                *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
    2650 
     2837       
    26512838        }
    26522839}
    … …  
    26582845
    26592846        /*
    2660          * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    2661          * For grid i, Lpi can be expressed in the basic coordinate system
    2662          * by:
    2663          *       Lpi_basic=[ h    0    0   0]
    2664          *               [ 0    h    0   0]
    2665          *               [ h    0    0   0]
    2666          *               [ 0    h    0   0]
    2667          *               [ 0    0    h   0]
    2668          *               [ 0    0    h   0]
    2669          *               [ 0    0  dh/dz 0]
    2670          *               [ 0    0  dh/dz 0]
    2671          *               [ 0    0  dh/dz 0]
    2672          *               [dh/dz 0  dh/dx 0]
    2673          *               [ 0 dh/dz dh/dy 0]
    2674          *               [ 0    0    0   h]
    2675          *               [ 0    0    0   h]
    2676          *               [ 0    0    0   h]
    2677          * where h is the interpolation function for grid i.
    2678          */
     2847        * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     2848        * For grid i, Lpi can be expressed in the basic coordinate system
     2849        * by:
     2850        *       Lpi_basic=[ h    0    0   0]
     2851        *                [ 0    h    0   0]
     2852        *                [ h    0    0   0]
     2853        *                [ 0    h    0   0]
     2854        *                [ 0    0    h   0]
     2855        *                [ 0    0    h   0]
     2856        *                [ 0    0  dh/dz 0]
     2857        *                [ 0    0  dh/dz 0]
     2858        *                [ 0    0  dh/dz 0]
     2859        *                [dh/dz 0  dh/dx 0]
     2860        *                [ 0 dh/dz dh/dy 0]
     2861        *                [ 0    0    0   h]
     2862        *                [ 0    0    0   h]
     2863        *                [ 0    0    0   h]
     2864        * where h is the interpolation function for grid i.
     2865        */
    26792866
    26802867        int i;
    … …  
    26902877        l1l2l3[1]=gauss_coord_tria[1];
    26912878        l1l2l3[2]=gauss_coord_tria[2];
    2692 
     2879       
    26932880        GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
    26942881
    2695 #ifdef _DELUG_
     2882        #ifdef _DELUG_
    26962883        for (i=0;i<3;i++){
    26972884                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
    26982885        }
    2699 #endif
     2886        #endif
    27002887
    27012888        /*Build LprimeStokes: */
    … …  
    27572944                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
    27582945                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
    2759 
    2760         }
    2761 
     2946       
     2947        }
     2948       
    27622949}
    27632950
    … …  
    27652952#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasicStokes"
    27662953void Penta::GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord){
    2767 
     2954       
    27682955        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    27692956         * basic coordinate system: */
    … …  
    28012988#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParamsStokes"
    28022989void Penta::GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord){
    2803 
     2990       
    28042991        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    28052992         * natural coordinate system) at the gaussian point. */
    … …  
    28933080        int     area_order=5;
    28943081        int       num_vert_gauss=5;
    2895 
     3082       
    28963083        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    28973084        double  viscosity;
    … …  
    29163103        double     tBD[27][8];
    29173104        double     P_terms[numdof];
    2918 
     3105       
    29193106        ParameterInputs* inputs=NULL;
    29203107        Tria*            tria=NULL;
    … …  
    29413128
    29423129        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    2943                 get tria gaussian points as well as segment gaussian points. For tria gaussian
    2944                 points, order of integration is 2, because we need to integrate the product tB*D*B'
    2945                 which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    2946                 points, same deal, which yields 3 gaussian points.*/
     3130           get tria gaussian points as well as segment gaussian points. For tria gaussian
     3131           points, order of integration is 2, because we need to integrate the product tB*D*B'
     3132           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3133           points, same deal, which yields 3 gaussian points.*/
    29473134
    29483135        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    … …  
    29643151                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    29653152                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    2966 
     3153               
    29673154                        /* Get Jacobian determinant: */
    29683155                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    … …  
    29843171                        GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord);
    29853172                        GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord);
    2986 
     3173                       
    29873174                        /*Get bubble part of Bprime */
    29883175                        for(i=0;i<8;i++){
    … …  
    30333220                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
    30343221                        gauss_coord[3]=-1;
    3035 
     3222                       
    30363223                        gauss_coord_tria[0]=*(first_gauss_area_coord+igarea);
    30373224                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
    … …  
    30403227                        /*Get the Jacobian determinant */
    30413228                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
    3042 
     3229                       
    30433230                        /* Get bed at gaussian point */
    30443231                        GetParameterValue(&bed,&b[0],gauss_coord);
    … …  
    30463233                        /*Get L matrix : */
    30473234                        tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1);
    3048 
     3235                               
    30493236                        /*Get water_pressure at gaussian point */
    30503237                        water_pressure=gravity*rho_water*bed;
    … …  
    30523239                        /*Get normal vecyor to the bed */
    30533240                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
    3054 
     3241                       
    30553242                        bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
    30563243                        bed_normal[1]=-surface_normal[1];
    … …  
    30643251                }
    30653252        } //if ( (onbed==1) && (shelf==1))
    3066 
     3253       
    30673254        /*Reduce the matrix */
    30683255        ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]);
    … …  
    30803267#define __FUNCT__ "Penta::ReduceVectorStokes"
    30813268void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){
    3082 
     3269               
    30833270        int i,j;
    30843271
    … …  
    31233310#define __FUNCT__ "Penta::GetNodalFunctionsStokes"
    31243311void Penta::GetNodalFunctionsStokes(double* l1l7, double* gauss_coord){
    3125 
     3312       
    31263313        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    31273314
    … …  
    31813368        int     dofs[3]={0,1,2};
    31823369        double  dt;
    3183         int     dt_exists;
    31843370
    31853371        double  vxvyvz_list[numgrids][3];
    … …  
    31953381
    31963382        /*matrices: */
    3197         double  K_terms[numdof][numdof]={0.0};
    3198         double  Ke_gaussian_conduct[numdof][numdof];
    3199         double  Ke_gaussian_advec[numdof][numdof];
    3200         double  Ke_gaussian_transient[numdof][numdof];
    3201         double  B[3][numdof];
    3202         double  Bprime[3][numdof];
    3203         double  B_conduct[3][numdof];
    3204         double  B_advec[3][numdof];
    3205         double  Bprime_advec[3][numdof];
    3206         double  L[numdof];
    3207         double  D_scalar;
    3208         double  D[3][3];
    3209         double  l1l2l3[3];
    3210         double  tl1l2l3D[3];
    3211         double  tBD[3][numdof];
    3212         double  tBD_conduct[3][numdof];
    3213         double  tBD_advec[3][numdof];
    3214         double  tLD[numdof];
    3215 
    3216         double  Jdet;
     3383        double     K_terms[numdof][numdof]={0.0};
     3384        double     Ke_gaussian_conduct[numdof][numdof];
     3385        double     Ke_gaussian_advec[numdof][numdof];
     3386        double     Ke_gaussian_transient[numdof][numdof];
     3387        double     B[3][numdof];
     3388        double     Bprime[3][numdof];
     3389        double     B_conduct[3][numdof];
     3390        double     B_advec[3][numdof];
     3391        double     Bprime_advec[3][numdof];
     3392        double     L[numdof];
     3393        double     D_scalar;
     3394        double     D[3][3];
     3395        double     l1l2l3[3];
     3396        double     tl1l2l3D[3];
     3397        double     tBD[3][numdof];
     3398        double     tBD_conduct[3][numdof];
     3399        double     tBD_advec[3][numdof];
     3400        double     tLD[numdof];
     3401
     3402        double     Jdet;
    32173403
    32183404        /*Material properties: */
    3219         double  gravity,rho_ice,rho_water;
    3220         double  heatcapacity,thermalconductivity;
    3221         double  mixed_layer_capacity,thermal_exchange_velocity;
     3405        double         gravity,rho_ice,rho_water;
     3406        double         heatcapacity,thermalconductivity;
     3407        double         mixed_layer_capacity,thermal_exchange_velocity;
    32223408
    32233409        /*Collapsed formulation: */
    … …  
    32323418        GetDofList(&doflist[0],&numberofdofspernode);
    32333419       
    3234         /*recovre material parameters: */
     3420        // /*recovre material parameters: */
    32353421        rho_water=matpar->GetRhoWater();
    32363422        rho_ice=matpar->GetRhoIce();
    … …  
    32383424        heatcapacity=matpar->GetHeatCapacity();
    32393425        thermalconductivity=matpar->GetThermalConductivity();
     3426
    32403427               
    32413428        /*recover extra inputs from users, dt and velocity: */
    … …  
    32943481                        D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar;
    32953482
     3483
    32963484                        /*  Do the triple product B'*D*B: */
    32973485                        MatrixMultiply(&B_conduct[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_conduct[0][0],0);
    … …  
    33033491                        GetB_advec(&B_advec[0][0],&xyz_list[0][0],gauss_coord);
    33043492                        GetBprime_advec(&Bprime_advec[0][0],&xyz_list[0][0],gauss_coord);
     3493
    33053494
    33063495                        //Build the D matrix
    … …  
    33223511                        MatrixMultiply(&B_advec[0][0],3,numdof,1,&D[0][0],3,3,0,&tBD_advec[0][0],0);
    33233512                        MatrixMultiply(&tBD_advec[0][0],numdof,3,0,&Bprime_advec[0][0],3,numdof,0,&Ke_gaussian_advec[0][0],0);
     3513
    33243514
    33253515                        /*Transient: */
    … …  
    33613551        xfree((void**)&vert_gauss_weights);
    33623552        xfree((void**)&vert_gauss_coord);
    3363 
     3553       
    33643554        //Ice/ocean heat exchange flux on ice shelf base
    33653555        if(onbed && shelf){
    … …  
    33853575         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    33863576         */
    3387 
     3577       
    33883578        int i;
    33893579        int calculationdof=3;
    … …  
    34193609         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    34203610         */
    3421 
     3611       
    34223612        int i;
    34233613        int calculationdof=3;
    … …  
    34543644         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    34553645         */
    3456 
     3646       
    34573647        int i;
    34583648        int calculationdof=3;
    … …  
    34773667#define __FUNCT__ "Penta::CreateKMatrixMelting"
    34783668void  Penta::CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    3479 
     3669       
    34803670        Tria* tria=NULL;
    34813671        if (!onbed){
    … …  
    34833673        }
    34843674        else{
    3485 
     3675               
    34863676                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    34873677                tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type);
    … …  
    35083698        /*Grid data: */
    35093699        double        xyz_list[numgrids][3];
    3510 
     3700       
    35113701        /* gaussian points: */
    35123702        int     num_area_gauss,igarea,igvert;
    … …  
    35213711        int area_order=2;
    35223712        int     num_vert_gauss=3;
    3523 
     3713       
    35243714        double         dt;
    35253715        double         vx_list[numgrids];
    … …  
    35373727        /* element parameters: */
    35383728        int     friction_type;
    3539 
     3729       
    35403730        int            dofs[3]={0,1,2};
    35413731        int            dofs1[1]={0};
    … …  
    35723762        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    35733763        GetDofList(&doflist[0],&numberofdofspernode);
    3574 
     3764       
    35753765        // /*recovre material parameters: */
    35763766        rho_water=matpar->GetRhoWater();
    … …  
    35923782                if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
    35933783        }
    3594 
     3784       
    35953785        for(i=0;i<numgrids;i++){
    35963786                vx_list[i]=vxvyvz_list[i][0];
    … …  
    36003790
    36013791        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    3602                 get tria gaussian points as well as segment gaussian points. For tria gaussian
    3603                 points, order of integration is 2, because we need to integrate the product tB*D*B'
    3604                 which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    3605                 points, same deal, which yields 3 gaussian points.: */
    3606 
     3792           get tria gaussian points as well as segment gaussian points. For tria gaussian
     3793           points, order of integration is 2, because we need to integrate the product tB*D*B'
     3794           which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3795           points, same deal, which yields 3 gaussian points.: */
     3796       
    36073797        /*Get gaussian points: */
    36083798        GaussPenta( &num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights,&vert_gauss_coord, &vert_gauss_weights, area_order, num_vert_gauss);
    … …  
    36193809                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
    36203810                        gauss_coord[3]=*(vert_gauss_coord+igvert);
    3621 
     3811       
    36223812                        /*Compute strain rate and viscosity: */
    36233813                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    3624                         matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     3814                        matice->GetViscosity3d(&viscosity,&epsilon[0]);
    36253815
    36263816                        /* Get Jacobian determinant: */
    … …  
    36613851        /* Ice/ocean heat exchange flux on ice shelf base */
    36623852        if(onbed && shelf){
     3853
    36633854
    36643855                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.

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

    r483 r586  
    110110                void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4);
    111111                void  GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4);
     112                void  VelocityExtrudeAllElements(Vec ug,double* ug_serial);
    112113                void  VelocityExtrude(Vec ug,double* ug_serial);
     114                void  ThicknessExtrude(Vec ug,double* ug_serial);
     115                void  VelocityDepthAverageAtBase(Vec ug,double* ug_serial);
    113116                void  SlopeExtrude(Vec sg,double* sg_serial);
    114117                void  ComputePressure(Vec p_g);

  • issm/trunk/src/c/objects/Tria.cpp

    r511 r586  
    1717#include "../include/typedefs.h"
    1818
     19
    1920/*For debugging purposes: */
    2021#define ELID 141 //element for which to print debug statements
    … …  
    3334Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],double tria_melting[3],
    3435                                double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_meanvel,double tria_epsvel,
    35                                 double tria_viscosity_overshoot){
     36                                double tria_viscosity_overshoot,int tria_artdiff){
    3637       
    3738        int i;
    … …  
    5051                melting[i]=tria_melting[i];
    5152                accumulation[i]=tria_accumulation[i];
    52                 geothermalflux[i]=tria_geothermalflux[i];
    53         }
     53                geothermalflux[i]=tria_geothermalflux[i]; }
    5454        matice=NULL;
    5555        matice_offset=UNDEF;
    … …  
    6464        onbed=1;
    6565        viscosity_overshoot=tria_viscosity_overshoot;
     66        artdiff=tria_artdiff;
    6667       
    6768        return;
    … …  
    9697        printf("   onbed: %i\n",onbed);
    9798        printf("   viscosity_overshoot=%g\n",viscosity_overshoot);
     99        printf("   artdiff=%g\n",artdiff);
    98100        printf("   nodes: \n");
    99101        if(nodes[0])nodes[0]->Echo();
    … …  
    146148        memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
    147149        memcpy(marshalled_dataset,&viscosity_overshoot,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
     150        memcpy(marshalled_dataset,&artdiff,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff);
    148151       
    149152        *pmarshalled_dataset=marshalled_dataset;
    … …  
    177180                +sizeof(epsvel)
    178181                +sizeof(viscosity_overshoot)
     182                +sizeof(artdiff)
    179183                +sizeof(int); //sizeof(int) for enum type
    180184}
    … …  
    220224        memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
    221225        memcpy(&viscosity_overshoot,marshalled_dataset,sizeof(viscosity_overshoot));marshalled_dataset+=sizeof(viscosity_overshoot);
     226        memcpy(&artdiff,marshalled_dataset,sizeof(artdiff));marshalled_dataset+=sizeof(artdiff);
    222227
    223228        /*nodes and materials are not pointing to correct objects anymore:*/
    … …  
    242247}
    243248
     249
    244250#undef __FUNCT__
    245251#define __FUNCT__ "Tria::Configure"
    … …  
    290296
    291297        }
     298        else if (analysis_type==PrognosticAnalysisEnum()){
     299
     300                CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
     301
     302        }
    292303        else{
    293304                throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
    294305        }
    295 }
     306
     307}
     308
    296309
    297310#undef __FUNCT__
    … …  
    299312
    300313void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     314
    301315
    302316        /* local declarations */
    … …  
    509523
    510524}
     525#undef __FUNCT__
     526#define __FUNCT__ "Tria::CreateKMatrixPrognostic"
     527void  Tria::CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     528
     529
     530        /* local declarations */
     531        int             i,j;
     532
     533        /* node data: */
     534        const int    numgrids=3;
     535        const int    NDOF1=1;
     536        const int    numdof=NDOF1*numgrids;
     537        double       xyz_list[numgrids][3];
     538        int          doflist[numdof];
     539        int          numberofdofspernode;
     540       
     541        /* gaussian points: */
     542        int     num_gauss,ig;
     543        double* first_gauss_area_coord  =  NULL;
     544        double* second_gauss_area_coord =  NULL;
     545        double* third_gauss_area_coord  =  NULL;
     546        double* gauss_weights           =  NULL;
     547        double  gauss_weight;
     548        double  gauss_l1l2l3[3];
     549
     550        /* matrices: */
     551        double L[numgrids];
     552        double DL[2][2]={0.0};
     553        double DLprime[2][2]={0.0};
     554        double DL_scalar;
     555        double Ke_gg[numdof][numdof]={0.0};//local element stiffness matrix
     556        double Ke_gg_gaussian[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     557        double Ke_gg_thickness[numdof][numdof]={0.0}; //stiffness matrix evaluated at the gaussian point.
     558       
     559        double Jdettria;
     560       
     561        /*input parameters for structural analysis (diagnostic): */
     562        double  vxvy_list[numgrids][2]={0,0};
     563        double  vx_list[numgrids]={0,0};
     564        double  vy_list[numgrids]={0,0};
     565        double  dvx[2];
     566        double  dvy[2];
     567        double  vx,vy;
     568        double  v_gauss[2]={0.0};
     569        double  K[2][2]={0.0};
     570        double  dt;
     571        int     dofs[2]={0,1};
     572        int     found=0;
     573
     574        ParameterInputs* inputs=NULL;
     575
     576        /*recover pointers: */
     577        inputs=(ParameterInputs*)vinputs;
     578
     579        /*recover extra inputs from users, at current convergence iteration: */
     580        found=inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     581        if(!found)throw ErrorException(__FUNCT__," could not find velocity_average  in inputs!");
     582
     583        for(i=0;i<numgrids;i++){
     584                vx_list[i]=vxvy_list[i][0];
     585                vy_list[i]=vxvy_list[i][1];
     586        }
     587       
     588        found=inputs->Recover("dt",&dt);
     589        if(!found)throw ErrorException(__FUNCT__," could not find dt in inputs!");
     590
     591        /* Get node coordinates and dof list: */
     592        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     593        GetDofList(&doflist[0],&numberofdofspernode);
     594
     595        //Create Artificial diffusivity once for all if requested
     596        if(artdiff){
     597                //Get the Jacobian determinant
     598                gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0;
     599                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     600
     601                //Build K matrix (artificial diffusivity matrix)
     602                v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]);
     603                v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]);
     604               
     605                K[0][0]=pow(Jdettria,.5)/2.0*fabs(v_gauss[0]);
     606                K[1][1]=pow(Jdettria,.5)/2.0*fabs(v_gauss[1]);
     607        }
     608
     609        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     610        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     611
     612        /* Start  looping on the number of gaussian points: */
     613        for (ig=0; ig<num_gauss; ig++){
     614                /*Pick up the gaussian point: */
     615                gauss_weight=*(gauss_weights+ig);
     616                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     617                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     618                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     619
     620                /* Get Jacobian determinant: */
     621                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3);
     622
     623                /*Get L matrix: */
     624                GetL(&L[0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     625
     626
     627                DL_scalar=gauss_weight*Jdettria;
     628
     629                /*  Do the triple product tL*D*L: */
     630                TripleMultiply( &L[0],1,numdof,1,
     631                                          &DL_scalar,1,1,0,
     632                                          &L[0],1,numdof,0,
     633                                          &Ke_gg_gaussian[0][0],0);
     634               
     635               
     636                /*Get B  and B prime matrix: */
     637                //GetB_prog(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
     638                //GetBprime_prog(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);
     639               
     640                //Get vx, vy and their derivatives at gauss point
     641                GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
     642                GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
     643               
     644                GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
     645                GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
     646
     647                //dvxdx=dvx[0];
     648                //dvydy=dvy[1];
     649
     650                DL_scalar=dt*gauss_weight*Jdettria;
     651
     652                //Create DL and DLprime matrix
     653                //DL[0][0]=DL_scalar*dvxdx;
     654                //DL[1][1]=DL_scalar*dvydy;
     655
     656                DLprime[0][0]=DL_scalar*vx;
     657                DLprime[1][1]=DL_scalar*vy;
     658
     659                //Do the triple product tL*D*L.
     660                //Ke_gg_thickness=B'*DL*B+B'*DLprime*Bprime;
     661
     662                /*TripleMultiply( &B[0][0],numdof,numdof,1,
     663                                          &DL_scalar,1,1,0,
     664                                          &L[0],1,numdof,0,
     665                                          &Ke_gg_gaussian[0][0],0);*/
     666
     667
     668
     669
     670
     671                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     672                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     673               
     674                #ifdef _DEBUGELEMENTS_
     675                if(my_rank==RANK && id==ELID){
     676                        printf("      B:\n");
     677                        for(i=0;i<3;i++){
     678                                for(j=0;j<numdof;j++){
     679                                        printf("%g ",B[i][j]);
     680                                }
     681                                printf("\n");
     682                        }
     683                        printf("      Bprime:\n");
     684                        for(i=0;i<3;i++){
     685                                for(j=0;j<numdof;j++){
     686                                        printf("%g ",Bprime[i][j]);
     687                                }
     688                                printf("\n");
     689                        }
     690                }
     691                #endif
     692        } // for (ig=0; ig<num_gauss; ig++)
     693
     694        /*Add Ke_gg to global matrix Kgg: */
     695        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     696
     697
     698        /*Do not forget to include friction: */
     699        if(!shelf){
     700                //CreateKMatrixPrognosticFriction(Kgg,inputs,analysis_type,sub_analysis_type);
     701        }
     702
     703
     704        #ifdef _DEBUGELEMENTS_
     705        if(my_rank==RANK && id==ELID){
     706                printf("      Ke_gg erms:\n");
     707                for( i=0; i<numdof; i++){
     708                        for (j=0;j<numdof;j++){
     709                                printf("%g ",Ke_gg[i][j]);
     710                        }
     711                        printf("\n");
     712                }
     713                printf("      Ke_gg row_indices:\n");
     714                for( i=0; i<numdof; i++){
     715                        printf("%i ",doflist[i]);
     716                }
     717
     718        }
     719        #endif
     720
     721        cleanup_and_return:
     722        xfree((void**)&first_gauss_area_coord);
     723        xfree((void**)&second_gauss_area_coord);
     724        xfree((void**)&third_gauss_area_coord);
     725        xfree((void**)&gauss_weights);
     726
     727}
     728
    511729
    512730#undef __FUNCT__

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

    r483 r586  
    4949                int    onbed;
    5050                double viscosity_overshoot;
     51                int    artdiff;
    5152
    5253        public:
    … …  
    5455                Tria();
    5556                Tria(int id,int mid,int mparid,int node_ids[3],double h[3],double s[3],double b[3],double k[3],double melting[3],double accumulation[3],double geothermalflux[3],
    56                                 int friction_type,double p,double q,int shelf,double meanvel,double epsvel,double viscosity_overshoot);
     57                                int friction_type,double p,double q,int shelf,double meanvel,double epsvel,double viscosity_overshoot,int artdiff);
    5758                ~Tria();
    5859
    … …  
    113114                void  CreatePVectorThermalShelf( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
    114115                void  CreatePVectorThermalSheet( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type);
     116                void  CreateKMatrixPrognostic(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type);
    115117
    116118

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

    r387 r586  
    3737#include "./OptPars.h"
    3838
     39
    3940#endif

  • issm/trunk/src/c/parallel/CreateFemModel.cpp

    r557 r586  
    6363        _printf_("   normalizing rigid body constraints matrix:\n");
    6464        NormalizeConstraintsx(&Gmn, Rmg,nodesets);
    65        
     65
    6666        _printf_("   configuring element and loads:\n");
    6767        ConfigureObjectsx(elements, loads, nodes, materials);
    68        
     68
    6969        _printf_("   process parameters:\n");
    7070        ProcessParamsx( parameters, partition);
    71        
     71
    7272        _printf_("   free ressources:\n");
    7373        DeleteModel(&model);
     74
    7475
    7576        /*Assign output pointers:*/

  • issm/trunk/src/c/parallel/diagnostic.cpp

    r472 r586  
    2323        char* outputfilename=NULL;
    2424        char* lockname=NULL;
     25        char* qmuname=NULL;
    2526        int   numberofnodes;
     27        int   qmu_analysis=0;
    2628
    2729        /*Fem models : */
    … …  
    5254        outputfilename=argv[3];
    5355        lockname=argv[4];
     56        qmuname=argv[5];
    5457
    5558        /*Open handle to data on disk: */
    … …  
    6871        CreateFemModel(&femmodels[4],fid,"slope_compute","");
    6972
     73
    7074        _printf_("initialize inputs:\n");
    7175        femmodels[0].parameters->FindParam((void*)&u_g_initial,"u_g");
    … …  
    7983        femmodels[0].parameters->DeleteObject((Object*)param);
    8084
    81         _printf_("call computational core:\n");
    82         diagnostic_core(&u_g,&p_g,femmodels,inputs);
     85        /*are we running the solutoin sequence, or a qmu wrapper around it? : */
     86        femmodels[0].parameters->FindParam((void*)&qmu_analysis,"qmu_analysis");
     87        if(!qmu_analysis){
     88                /*run diagnostic analysis: */
     89                _printf_("call computational core:\n");
     90                diagnostic_core(&u_g,&p_g,femmodels,inputs);
    8391
    84         _printf_("write results to disk:\n");
    85         OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename);
     92                _printf_("write results to disk:\n");
     93                OutputDiagnostic(u_g,p_g,&femmodels[0],outputfilename);
     94        }
     95        else{
     96       
     97                /*run qmu analysis: */
     98                _printf_("calling qmu analysis on diagnostic core:\n");
     99               
     100                qmu(qmuname,&femmodels[0],inputs,DiagnosticAnalysisEnum(),NoneAnalysisEnum());
     101        }
    86102
    87103        _printf_("write lock file:\n");

  • issm/trunk/src/c/parallel/diagnostic_core.cpp

    r465 r586  
    1414
    1515void diagnostic_core(Vec* pug, Vec* ppg,FemModel* fems, ParameterInputs* inputs){
     16
     17        extern int my_rank;
    1618
    1719        /*fem models: */

  • issm/trunk/src/c/parallel/diagnostic_core_nonlinear.cpp

    r551 r586  
    4646        kflag=1; pflag=1;
    4747
     48
     49
    4850        fem->parameters->FindParam((void*)&connectivity,"connectivity");
    4951        fem->parameters->FindParam((void*)&numberofdofspernode,"numberofdofspernode");
    … …  
    5759        /*Copy loads for backup: */
    5860        loads=fem->loads->Copy();
    59        
     61
    6062        /*Initialize ug and ug_old */
    6163        if (numberofdofspernode>=3)dofs[2]=1;//only keep vz if running with more than 3 dofs per node

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

    r568 r586  
    3131void OutputThermal(Vec* t_g,Vec* m_g, double* time,FemModel* femmodels,char* filename);
    3232void OutputControl(Vec u_g,double* p_g, double* J, int nsteps, Vec partition,char* filename,NodeSets* nodesets);
     33void OutputPrognostic(Vec h_g,FemModel* femmodel,char* filename);
    3334void WriteLockFile(char* filename);
    3435
    3536void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,char* analysis_type,char* sub_analysis_type);
    3637//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femmodel,char* filename);
     38void qmu(const char* dakota_input_file,FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     39void SpawnCore(double* responses,double* variables,char** variable_descriptors,int numvariables, FemModel* femmodels,ParameterInputs* inputs,int analysis_type,int sub_analysis_type);
     40void DakotaResponses(double* responses,char** responses_descriptors,int numresponses,Vec u_g,Vec p_g,int analysis_type,int sub_analysis_type);
    3741
    3842#endif

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

    r496 r586  
    4242                numdofs=1;
    4343        }
     44        else if (analysis_type==PrognosticAnalysisEnum()){
     45                numdofs=1;
     46        }
    4447        else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis type: ",analysis_type,"  not implemented yet!"));
    4548

  • issm/trunk/src/c/toolkits/metis/metisincludes.h

    r1 r586  
    99#include <metis.h>
    1010
     11extern "C" void METIS_PartMeshNodal(int *, int *, idxtype *, int *, int *, int *, int *, idxtype *, idxtype *);
     12
    1113#endif

  • issm/trunk/src/c/toolkits/petsc/patches/petscpatches.h

    r434 r586  
    4242void MatToSerial(double** poutmatrix,Mat matrix);
    4343void VecDuplicatePatch(Vec* output, Vec input);
     44Vec  SerialToVec(double* vector,int vector_size);
    4445
    4546

  • issm/trunk/src/m/classes/@model/model.m

    r465 r586  
    274274        md.dakotaout='';
    275275        md.dakotadat='';
     276        md.qmu_analysis=0;
     277        md.iresp=1;
    276278
    277279        md.npart=0;

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

    r465 r586  
    2727end
    2828
    29 fprintf(fid,' %s %s.bin %s.outbin %s.lock 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name);
     29fprintf(fid,' %s %s.bin %s.outbin %s.lock qmu.in 2> %s.errlog >%s.outlog & ',executionpath,md.name,md.name,md.name,md.name,md.name);
    3030fclose(fid);

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

    r1 r586  
    2020disp('uploading input file and queueing script');
    2121if strcmpi(hostname,md.cluster),
    22         system(['cp ' md.name '.bin' ' ' md.name '.queue' ' ' executionpath]);
     22        if md.qmu_analysis,
     23                system(['cp ' md.name '.bin' ' ' md.name '.queue qmu/qmu.in ' ' ' executionpath]);
     24        else
     25                system(['cp ' md.name '.bin' ' ' md.name '.queue' ' ' executionpath]);
     26        end
    2327else
    24         system(['scp ' md.name '.bin' ' ' md.name '.queue' ' ' md.cluster ':' executionpath]);
     28        if md.qmu_analysis,
     29                system(['scp ' md.name '.bin' ' ' md.name '.queue qmu/qmu.in ' ' ' md.cluster ':' executionpath]);
     30        else
     31                system(['scp ' md.name '.bin' ' ' md.name '.queue' ' ' md.cluster ':' executionpath]);
     32        end
    2533end
    2634

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

    r574 r586  
    153153        WriteData(fid,md.rifts(i).friction,'Scalar',['friction' num2str(i)]);
    154154end
     155
     156%Qmu
     157WriteData(fid,md.qmu_analysis,'Integer','qmu_analysis');
     158if md.qmu_analysis,
     159        responses=md.responses(md.iresp).rf;
     160        WriteData(fid,numel(responses),'Integer','numberofresponses');
     161        for i=1:numel(responses),
     162                WriteData(fid,responses(i).descriptor,'String',['descriptor' num2str(i-1)]);
     163        end
     164end
    155165       
    156166%Get penalties to connect collapsed and non-collapsed 3d meshes:

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

    r490 r586  
    4646end
    4747if ~found,
    48         sub_analysis_type='';
     48        sub_analysis_type='none';
    4949end
    5050
    … …  
    5656        strcmpi(analysis_type,'mesh') |  ...
    5757        strcmpi(analysis_type,'mesh2grid') |  ...
    58         strcmpi(analysis_type,'qmu') |  ...
    5958        strcmpi(analysis_type,'transient') ),
    6059        error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']);
    6160end
    62 if ~(strcmpi(sub_analysis_type,'') |  ...
     61if ~(strcmpi(sub_analysis_type,'none') |  ...
    6362        strcmpi(sub_analysis_type,'steady') |  ...
    6463        strcmpi(sub_analysis_type,'horiz') |  ...

  • issm/trunk/src/m/solutions/dakota/dakota_m_write.m

    r449 r586  
    208208
    209209disp(sprintf('  Writing %d %s responses.',length(dresp),class(dresp)));
    210 
    211210for i=1:length(dresp)
    212211    ifnvals=ifnvals+1;

  • issm/trunk/src/m/utils/Nightly/testsgetpackage.m

    r522 r586  
    2525elseif strcmpi(string,'cielo_parallel'),
    2626        package='cielo';
    27         md.cluster='wilkes';
     27        md.cluster='astrid';
    2828
    2929else

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

    r358 r586  
    3838                                SystemMatrices\
    3939                                Test\
     40                                ThicknessExtrude\
    4041                                TriMesh\
    4142                                TriMeshProcessRifts\
    4243                                TriMeshRefine\
    4344                                UpdateFromInputs\
    44                                 VelocityExtrude
     45                                VelocityExtrude\
     46                                VelocityDepthAverage
    4547
    4648endif
    … …  
    167169                          SystemMatrices/SystemMatrices.h
    168170
     171ThicknessExtrude_SOURCES = ThicknessExtrude/ThicknessExtrude.cpp\
     172                          ThicknessExtrude/ThicknessExtrude.h
     173
    169174TriMesh_SOURCES = TriMesh/TriMesh.cpp\
    170175                          TriMesh/TriMesh.h
    … …  
    182187                          VelocityExtrude/VelocityExtrude.h
    183188
     189VelocityDepthAverage_SOURCES = VelocityDepthAverage/VelocityDepthAverage.cpp\
     190                          VelocityDepthAverage/VelocityDepthAverage.h
Note: See TracChangeset for help on using the changeset viewer.