Changeset 2711


Ignore:
Timestamp:
12/09/09 07:37:38 (15 years ago)
Author:
Mathieu Morlighem
Message:

Added folds in Penta

Location:
issm/trunk/src/c/objects
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r2671 r2711  
    1616#include "../include/typedefs.h"
    1717
     18/*FUNCTION Penta constructor {{{1*/
    1819Penta::Penta(){
    1920        return;
    2021}
     22/*}}}*/
     23/*FUNCTION Penta creation {{{1*/
    2124Penta::Penta( int penta_id, int penta_mid, int penta_mparid, int penta_numparid, 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,  int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6],
    23                                  int penta_thermal_steadystate,bool penta_onwater){
    24        
     25                        double penta_p, double penta_q, int penta_shelf, int penta_onbed, int penta_onsurface,  int penta_collapse, double penta_melting[6], double penta_accumulation[6], double penta_geothermalflux[6],
     26                        int penta_thermal_steadystate,bool penta_onwater){
     27
    2528        int i;
    2629
    … …  
    6164        return;
    6265}
    63 
     66/*}}}*/
     67/*FUNCTION Penta destructor {{{1*/
    6468Penta::~Penta(){
    6569        return;
    6670}
     71/*}}}*/
     72/*FUNCTION Penta Echo {{{1*/
    6773void Penta::Echo(void){
    6874
    … …  
    8288        printf("   b=[%g,%g,%g,%g,%g,%g]\n",b[0],b[1],b[2],b[3],b[4],b[5]);
    8389        printf("   k=[%g,%g,%g,%g,%g,%g]\n",k[0],k[1],k[2],k[3],k[4],k[5]);
    84        
     90
    8591        printf("   friction_type: %i\n",friction_type);
    8692        printf("   p: %g\n",p);
    … …  
    9197        printf("   onsurface: %i\n",onsurface);
    9298        printf("   collapse: %i\n",collapse);
    93        
     99
    94100        printf("   melting=[%g,%g,%g,%g,%g,%g]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]);
    95101        printf("   accumulation=[%g,%g,%g,%g,%g,%g]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]);
    … …  
    98104        return;
    99105}
     106/*}}}*/
     107/*FUNCTION Penta DeepEcho {{{1*/
    100108void Penta::DeepEcho(void){
    101109
    … …  
    115123        printf("   b=[%i,%i,%i,%i,%i,%i]\n",b[0],b[1],b[2],b[3],b[4],b[5]);
    116124        printf("   k=[%i,%i,%i,%i,%i,%i]\n",k[0],k[1],k[2],k[3],k[4],k[5]);
    117        
     125
    118126        printf("   friction_type: %i\n",friction_type);
    119127        printf("   p: %g\n",p);
    … …  
    124132        printf("   onsurface: %i\n",onsurface);
    125133        printf("   collapse: %i\n",collapse);
    126        
     134
    127135        printf("   melting=[%i,%i,%i,%i,%i,%i]\n",melting[0],melting[1],melting[2],melting[3],melting[4],melting[5]);
    128136        printf("   accumulation=[%i,%i,%i,%i,%i,%i]\n",accumulation[0],accumulation[1],accumulation[2],accumulation[3],accumulation[4],accumulation[5]);
    … …  
    131139        return;
    132140}
     141/*}}}*/
     142/*FUNCTION Penta Marshall {{{1*/
    133143void  Penta::Marshall(char** pmarshalled_dataset){
    134144
    … …  
    141151        /*get enum type of Penta: */
    142152        enum_type=PentaEnum();
    143        
     153
    144154        /*marshall enum: */
    145155        memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);
    146        
     156
    147157        /*marshall Penta data: */
    148158        memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
    … …  
    175185        memcpy(marshalled_dataset,&geothermalflux,sizeof(geothermalflux));marshalled_dataset+=sizeof(geothermalflux);
    176186        memcpy(marshalled_dataset,&thermal_steadystate,sizeof(thermal_steadystate));marshalled_dataset+=sizeof(thermal_steadystate);
    177        
     187
    178188        *pmarshalled_dataset=marshalled_dataset;
    179189        return;
    180190}
    181                
     191/*}}}*/
     192/*FUNCTION Penta MarshallSize {{{1*/
    182193int   Penta::MarshallSize(){
    183194
    184195        return sizeof(id)+
    185                 sizeof(mid)+
    186                 sizeof(mparid)+
    187                 sizeof(node_ids)+
    188                 sizeof(nodes)+
    189                 sizeof(node_offsets)+
    190                 sizeof(matice)+
    191                 sizeof(matice_offset)+
    192                 sizeof(matpar)+
    193                 sizeof(matpar_offset)+
    194                 +sizeof(numparid)+
    195                 +sizeof(numpar)+
    196                 +sizeof(numpar_offset)+
    197                 sizeof(h)+
    198                 sizeof(s)+
    199                 sizeof(b)+
    200                 sizeof(k)+
    201                 sizeof(friction_type)+
    202                 sizeof(p)+
    203                 sizeof(q)+
    204                 sizeof(shelf)+
    205                 sizeof(onbed)+
    206                 sizeof(onwater)+
    207                 sizeof(onsurface)+
    208                 sizeof(collapse)+
    209                 sizeof(melting)+
    210                 sizeof(accumulation)+
    211                 sizeof(geothermalflux)+
    212                 sizeof(thermal_steadystate) +
    213                 sizeof(int); //sizeof(int) for enum type
    214 }
    215 
     196          sizeof(mid)+
     197          sizeof(mparid)+
     198          sizeof(node_ids)+
     199          sizeof(nodes)+
     200          sizeof(node_offsets)+
     201          sizeof(matice)+
     202          sizeof(matice_offset)+
     203          sizeof(matpar)+
     204          sizeof(matpar_offset)+
     205          +sizeof(numparid)+
     206          +sizeof(numpar)+
     207          +sizeof(numpar_offset)+
     208          sizeof(h)+
     209          sizeof(s)+
     210          sizeof(b)+
     211          sizeof(k)+
     212          sizeof(friction_type)+
     213          sizeof(p)+
     214          sizeof(q)+
     215          sizeof(shelf)+
     216          sizeof(onbed)+
     217          sizeof(onwater)+
     218          sizeof(onsurface)+
     219          sizeof(collapse)+
     220          sizeof(melting)+
     221          sizeof(accumulation)+
     222          sizeof(geothermalflux)+
     223          sizeof(thermal_steadystate) +
     224          sizeof(int); //sizeof(int) for enum type
     225}
     226/*}}}*/
     227/*FUNCTION Penta GetName {{{1*/
    216228char* Penta::GetName(void){
    217229        return "penta";
    218230}
    219 
     231/*}}}*/
     232/*FUNCTION Penta Demarshall {{{1*/
    220233void  Penta::Demarshall(char** pmarshalled_dataset){
    221234
    … …  
    269282        return;
    270283}
     284/*}}}*/
     285/*FUNCTION Penta Enum {{{1*/
    271286int Penta::Enum(void){
    272287
    … …  
    274289
    275290}
    276 int    Penta::GetId(void){ return id; }
    277 
     291/*}}}*/
     292/*FUNCTION Penta GetId {{{1*/
     293int    Penta::GetId(void){
     294        return id;
     295}
     296/*}}}*/
     297/*FUNCTION Penta MyRank {{{1*/
    278298int    Penta::MyRank(void){
    279299        extern int my_rank;
    280300        return my_rank;
    281301}
    282 
    283 
    284 #undef __FUNCT__
     302/*}}}*/
     303/*FUNCTION Penta Configure {{{1*/
     304#undef __FUConfigure NCT__
    285305#define __FUNCT__ "Penta::Configure"
    286306void  Penta::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin,void* pparametersin){
    287307
    288308        int i;
    289        
     309
    290310        DataSet* loadsin=NULL;
    291311        DataSet* nodesin=NULL;
    … …  
    301321        /*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
    302322        ResolvePointers((Object**)nodes,node_ids,node_offsets,6,nodesin);
    303        
     323
    304324        /*Same for materials: */
    305325        ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
    … …  
    310330
    311331}
    312 
     332/*}}}*/
     333/*FUNCTION Penta CreateKMatrix {{{1*/
    313334#undef __FUNCT__
    314335#define __FUNCT__ "Penta::CreateKMatrix"
    … …  
    323344        }
    324345        else if (analysis_type==DiagnosticAnalysisEnum()){
    325        
     346
    326347                if (sub_analysis_type==HorizAnalysisEnum()){
    327                
     348
    328349                        CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type,sub_analysis_type);
    329350                }
    330351                else if (sub_analysis_type==VertAnalysisEnum()){
    331                
     352
    332353                        CreateKMatrixDiagnosticVert( Kgg,inputs,analysis_type,sub_analysis_type);
    333354                }
    334355                else if (sub_analysis_type==StokesAnalysisEnum()){
    335                
     356
    336357                        CreateKMatrixDiagnosticStokes( Kgg,inputs,analysis_type,sub_analysis_type);
    337                
     358
    338359                }
    339360                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
    340361        }
    341362        else if (analysis_type==SlopeComputeAnalysisEnum()){
    342                
     363
    343364                CreateKMatrixSlopeCompute( Kgg,inputs,analysis_type,sub_analysis_type);
    344365        }
    345366        else if (analysis_type==PrognosticAnalysisEnum()){
    346                
     367
    347368                CreateKMatrixPrognostic( Kgg,inputs,analysis_type,sub_analysis_type);
    348369        }
    349370        else if (analysis_type==ThermalAnalysisEnum()){
    350                
     371
    351372                CreateKMatrixThermal( Kgg,inputs,analysis_type,sub_analysis_type);
    352373        }
    353374        else if (analysis_type==MeltingAnalysisEnum()){
    354                
     375
    355376                CreateKMatrixMelting( Kgg,inputs,analysis_type,sub_analysis_type);
    356377        }
    … …  
    360381
    361382}
    362 
     383/*}}}*/
     384/*FUNCTION Penta CreateKMatrixDiagnosticHoriz {{{1*/
    363385#undef __FUNCT__
    364386#define __FUNCT__ "Penta:CreateKMatrixDiagnosticHoriz"
    … …  
    375397        int          doflist[numdof];
    376398        int          numberofdofspernode;
    377        
    378        
     399
     400
    379401        /* 3d gaussian points: */
    380402        int     num_gauss,ig;
    … …  
    392414        int     num_area_gauss;
    393415        double  gauss_weight;
    394        
     416
    395417        /* 2d gaussian point: */
    396418        int     num_gauss2d;
    … …  
    424446        double Ke_gg_drag_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    425447        double Jdet;
    426        
     448
    427449        /*slope: */
    428450        double  slope[2]={0.0,0.0};
    … …  
    491513                }
    492514
    493                 #ifdef _DEBUGELEMENTS_
     515#ifdef _DEBUGELEMENTS_
    494516                if(my_rank==RANK && id==ELID){
    495517                        printf("El id %i Rank %i PentaElement input list before gaussian loop: \n",ELID,RANK);
    … …  
    508530                        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]);
    509531                }
    510                 #endif
     532#endif
    511533
    512534                /*Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    … …  
    521543                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);
    522544
    523                 #ifdef _DEBUGGAUSS_
     545#ifdef _DEBUGGAUSS_
    524546                if(my_rank==RANK && id==ELID){
    525547                        printf("El id %i Rank %i PentaElement gauss points\n",ELID,RANK);
    … …  
    531553                        }       
    532554                }
    533                 #endif
     555#endif
    534556                /* Start  looping on the number of gaussian points: */
    535557                for (ig1=0; ig1<num_area_gauss; ig1++){
    536558                        for (ig2=0; ig2<num_vert_gauss; ig2++){
    537                        
     559
    538560                                /*Pick up the gaussian point: */
    539561                                gauss_weight1=*(area_gauss_weights+ig1);
    540562                                gauss_weight2=*(vert_gauss_weights+ig2);
    541563                                gauss_weight=gauss_weight1*gauss_weight2;
    542                                
    543                                
     564
     565
    544566                                gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    545567                                gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    … …  
    551573                                GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
    552574                                GetStrainRate(&oldepsilon[0],&oldvxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3l4);
    553                        
     575
    554576                                /*Get viscosity: */
    555577                                matice->GetViscosity3d(&viscosity, &epsilon[0]);
    … …  
    562584                                /* Get Jacobian determinant: */
    563585                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
    564        
     586
    565587                                /*Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant
    566588                                  onto this scalar matrix, so that we win some computational time: */
    567                                
     589
    568590                                newviscosity=viscosity+numpar->viscosity_overshoot*(viscosity-oldviscosity);
    569591                                D_scalar=newviscosity*gauss_weight*Jdet;
    … …  
    571593                                        D[i][i]=D_scalar;
    572594                                }
    573                
     595
    574596                                /*  Do the triple product tB*D*Bprime: */
    575597                                TripleMultiply( &B[0][0],5,numdof,1,
    576                                                 &D[0][0],5,5,0,
    577                                                 &Bprime[0][0],5,numdof,0,
    578                                                 &Ke_gg_gaussian[0][0],0);
     598                                                        &D[0][0],5,5,0,
     599                                                        &Bprime[0][0],5,numdof,0,
     600                                                        &Ke_gg_gaussian[0][0],0);
    579601
    580602                                /* Add the Ke_gg_gaussian, and optionally Ke_gg_gaussian onto Ke_gg: */
    … …  
    586608                        } //for (ig2=0; ig2<num_vert_gauss; ig2++)
    587609                } //for (ig1=0; ig1<num_area_gauss; ig1++)
    588                
     610
    589611
    590612                /*Add Ke_gg to global matrix Kgg: */
    591613                MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    592        
     614
    593615                //Deal with 2d friction at the bedrock interface
    594616                if((onbed && !shelf)){
    … …  
    605627        }
    606628
    607         cleanup_and_return:
     629cleanup_and_return:
    608630        xfree((void**)&first_gauss_area_coord);
    609631        xfree((void**)&second_gauss_area_coord);
    … …  
    618640
    619641}
    620 
     642/*}}}*/
     643/*FUNCTION Penta CreateKMatrixDiagnosticVert {{{1*/
    621644#undef __FUNCT__
    622645#define __FUNCT__ "Penta:CreateKMatrixDiagnosticVert"
    … …  
    668691        /*recover pointers: */
    669692        inputs=(ParameterInputs*)vinputs;
    670        
     693
    671694
    672695        /*If this element  is on the surface, we have a dynamic boundary condition that applies, as a stiffness
    … …  
    701724
    702725        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);
    703         #ifdef _ISSM_DEBUG_
     726#ifdef _ISSM_DEBUG_
    704727        for (i=0;i<num_area_gauss;i++){
    705728                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));
    … …  
    708731                printf("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
    709732        }
    710         #endif
    711        
     733#endif
     734
    712735        /* Start  looping on the number of gaussian points: */
    713736        for (ig1=0; ig1<num_area_gauss; ig1++){
    714737                for (ig2=0; ig2<num_vert_gauss; ig2++){
    715                
     738
    716739                        /*Pick up the gaussian point: */
    717740                        gauss_weight1=*(area_gauss_weights+ig1);
    718741                        gauss_weight2=*(vert_gauss_weights+ig2);
    719742                        gauss_weight=gauss_weight1*gauss_weight2;
    720                        
     743
    721744                        gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    722745                        gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    … …  
    734757                        /*  Do the triple product tB*D*Bprime: */
    735758                        TripleMultiply( &B[0][0],1,numgrids,1,
    736                                         &DL_scalar,1,1,0,
    737                                         &Bprime[0][0],1,numgrids,0,
    738                                         &Ke_gg_gaussian[0][0],0);
     759                                                &DL_scalar,1,1,0,
     760                                                &Bprime[0][0],1,numgrids,0,
     761                                                &Ke_gg_gaussian[0][0],0);
    739762
    740763                        /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    … …  
    749772        /*Add Ke_gg to global matrix Kgg: */
    750773        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    751        
    752         cleanup_and_return:
     774
     775cleanup_and_return:
    753776        xfree((void**)&first_gauss_area_coord);
    754777        xfree((void**)&second_gauss_area_coord);
    … …  
    758781        xfree((void**)&vert_gauss_weights);
    759782}
    760 
     783/*}}}*/
     784/*FUNCTION Penta CreateKMatrixDiagnosticStokes {{{1*/
    761785#undef __FUNCT__
    762786#define __FUNCT__ "Penta:CreateKMatrixDiagnosticStokes"
    … …  
    786810        /*Grid data: */
    787811        double        xyz_list[numgrids][3];
    788        
     812
    789813        /*parameters: */
    790814        double             xyz_list_tria[numgrids2d][3];
    … …  
    810834        double     DLStokes[14][14]={0.0};
    811835        double     tLDStokes[numdof2d][14];
    812        
     836
    813837        /* gaussian points: */
    814838        int     num_area_gauss;
    … …  
    832856        double  alpha2_list[numgrids2d];
    833857        double  alpha2_gauss;
    834        
     858
    835859        ParameterInputs* inputs=NULL;
    836860
    … …  
    847871                }
    848872        }
    849        
     873
    850874        /*recovre material parameters: */
    851875        rho_water=matpar->GetRhoWater();
    … …  
    861885
    862886        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    863            get tria gaussian points as well as segment gaussian points. For tria gaussian
    864            points, order of integration is 2, because we need to integrate the product tB*D*B'
    865            which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    866            points, same deal, which yields 3 gaussian points.*/
     887                get tria gaussian points as well as segment gaussian points. For tria gaussian
     888                points, order of integration is 2, because we need to integrate the product tB*D*B'
     889                which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     890                points, same deal, which yields 3 gaussian points.*/
    867891
    868892        area_order=5;
    … …  
    888912                        /*Compute strain rate: */
    889913                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    890                
     914
    891915                        /*Get viscosity: */
    892916                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    … …  
    926950                /*Build alpha2_list used by drag stiffness matrix*/
    927951                Friction* friction=NewFriction();
    928                
     952
    929953                /*Initialize all fields: */
    930954                friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    931955                strcpy(friction->element_type,"2d");
    932                
     956
    933957                friction->gravity=matpar->GetG();
    934958                friction->rho_ice=matpar->GetRhoIce();
    … …  
    953977                        }
    954978                }
    955                
     979
    956980                xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
    957981                GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
    … …  
    964988                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
    965989                        gauss_coord[3]=-1;
    966                        
     990
    967991                        gauss_coord_tria[0]=*(first_gauss_area_coord+igarea);
    968992                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
    969993                        gauss_coord_tria[2]=*(third_gauss_area_coord+igarea);
    970        
     994
    971995                        /*Get the Jacobian determinant */
    972996                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
    973                        
     997
    974998                        /*Get L matrix if viscous basal drag present: */
    975999                        GetLStokes(&LStokes[0][0],  gauss_coord_tria);
    9761000                        GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss_coord_tria, gauss_coord);
    977                                
     1001
    9781002                        /*Compute strain rate: */
    9791003                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    980                
     1004
    9811005                        /*Get viscosity at last iteration: */
    9821006                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    983        
     1007
    9841008                        /*Get normal vecyor to the bed */
    9851009                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
    986                        
     1010
    9871011                        bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
    9881012                        bed_normal[1]=-surface_normal[1];
    … …  
    10181042                }
    10191043        } //if ( (onbed==1) && (shelf==0))
    1020        
     1044
    10211045        /*Reduce the matrix */
    10221046        ReduceMatrixStokes(&Ke_reduced[0][0], &Ke_temp[0][0]);
    … …  
    10301054        /*Add Ke_gg to global matrix Kgg: */
    10311055        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
    1032        
    1033         cleanup_and_return:
    1034        
     1056
     1057cleanup_and_return:
     1058
    10351059        /*Free ressources:*/
    10361060        xfree((void**)&first_gauss_area_coord);
    … …  
    10431067        return;
    10441068}
    1045 
     1069/*}}}*/
     1070/*FUNCTION Penta CreatePVector {{{1*/
    10461071#undef __FUNCT__
    10471072#define __FUNCT__ "Penta::CreatePVector"
    … …  
    10501075        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    10511076        if (analysis_type==ControlAnalysisEnum()){
    1052                
     1077
    10531078                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
    10541079        }
    10551080        else if (analysis_type==DiagnosticAnalysisEnum()){
    1056        
     1081
    10571082                if (sub_analysis_type==HorizAnalysisEnum()){
    10581083
    … …  
    10601085                }
    10611086                else if (sub_analysis_type==VertAnalysisEnum()){
    1062                
     1087
    10631088                        CreatePVectorDiagnosticVert( pg,inputs,analysis_type,sub_analysis_type);
    10641089                }
    10651090                else if (sub_analysis_type==StokesAnalysisEnum()){
    1066                
     1091
    10671092                        CreatePVectorDiagnosticStokes( pg,inputs,analysis_type,sub_analysis_type);
    10681093                }
    … …  
    10701095        }
    10711096        else if (analysis_type==SlopeComputeAnalysisEnum()){
    1072                
     1097
    10731098                CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
    10741099        }
    10751100        else if (analysis_type==PrognosticAnalysisEnum()){
    1076                
     1101
    10771102                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    10781103        }
    10791104        else if (analysis_type==ThermalAnalysisEnum()){
    1080                
     1105
    10811106                CreatePVectorThermal( pg,inputs,analysis_type,sub_analysis_type);
    10821107        }
    10831108        else if (analysis_type==MeltingAnalysisEnum()){
    1084                
     1109
    10851110                CreatePVectorMelting( pg,inputs,analysis_type,sub_analysis_type);
    10861111        }
    … …  
    10901115
    10911116}
     1117/*}}}*/
     1118/*FUNCTION Penta UpdateFromInputs {{{1*/
    10921119#undef __FUNCT__
    10931120#define __FUNCT__ "Penta::UpdateFromInputs"
    … …  
    11261153                }
    11271154        }
    1128        
     1155
    11291156        if(inputs->Recover("temperature_average",&temperature_list[0],1,dofs,6,(void**)nodes)){
    11301157                if(matice && collapse){
    … …  
    11361163                }
    11371164        }
    1138        
     1165
    11391166        if(inputs->Recover("B",&B_list[0],1,dofs,6,(void**)nodes)){
    11401167                if(matice){
    … …  
    11451172
    11461173}
    1147 
     1174/*}}}*/
     1175/*FUNCTION Penta GetMatPar {{{1*/
    11481176void* Penta::GetMatPar(){
    11491177        return matpar;
    11501178}
    1151 
     1179/*}}}*/
     1180/*FUNCTION Penta GetShelf {{{1*/
    11521181int   Penta::GetShelf(){
    11531182        return shelf;
    11541183}
    1155 
    1156 
     1184/*}}}*/
     1185/*FUNCTION Penta GetNodes {{{1*/
    11571186void  Penta::GetNodes(void** vpnodes){
    11581187        int i;
    … …  
    11631192        }
    11641193}
    1165                
     1194/*}}}*/
     1195/*FUNCTION Penta GetOnBed {{{1*/
    11661196int Penta::GetOnBed(){
    11671197        return onbed;
    11681198}
    1169 
    1170 void          Penta::GetThicknessList(double* thickness_list){
     1199/*}}}*/
     1200/*FUNCTION Penta GetThicknessList {{{1*/
     1201void Penta::GetThicknessList(double* thickness_list){
    11711202
    11721203        int i;
    11731204        for(i=0;i<6;i++)thickness_list[i]=h[i];
    11741205}
    1175 void          Penta::GetBedList(double* bed_list){
    1176        
     1206/*}}}*/
     1207/*FUNCTION Penta GetBedList {{{1*/
     1208void Penta::GetBedList(double* bed_list){
     1209
    11771210        int i;
    11781211        for(i=0;i<6;i++)bed_list[i]=b[i];
    11791212
    11801213}
    1181 
     1214/*}}}*/
     1215/*FUNCTION Penta copy {{{1*/
    11821216Object* Penta::copy() {
    11831217        return new Penta(*this);
    11841218}
    1185 
     1219/*}}}*/
     1220/*FUNCTION Penta Du {{{1*/
    11861221#undef __FUNCT__
    11871222#define __FUNCT__ "Penta::Du"
    … …  
    11931228        /*If on water, skip: */
    11941229        if(onwater)return;
    1195        
     1230
    11961231        /*Bail out if this element if:
    11971232         * -> Non collapsed and not on the surface
    … …  
    12101245        }
    12111246        else{
    1212                
     1247
    12131248                tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
    12141249                tria->Du(du_g,inputs,analysis_type,sub_analysis_type);
    … …  
    12171252        }
    12181253}
    1219 
     1254/*}}}*/
     1255/*FUNCTION Penta Gradj {{{1*/
    12201256#undef __FUNCT__
    12211257#define __FUNCT__ "Penta::Gradj"
    12221258void  Penta::Gradj(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type,char* control_type){
    1223        
     1259
    12241260        /*If on water, skip grad (=0): */
    12251261        if(onwater)return;
    … …  
    12331269        else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
    12341270}
    1235 
     1271/*}}}*/
     1272/*FUNCTION Penta GradjDrag {{{1*/
    12361273#undef __FUNCT__
    12371274#define __FUNCT__ "Penta::GradjDrag"
    12381275void  Penta::GradjDrag(Vec grad_g,void* inputs,int analysis_type,int sub_analysis_type){
    1239        
     1276
    12401277        Tria* tria=NULL;
    12411278
    … …  
    12671304        else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
    12681305}
    1269 
     1306/*}}}*/
     1307/*FUNCTION Penta GradjB {{{1*/
    12701308#undef __FUNCT__
    12711309#define __FUNCT__ "Penta::GradjB"
    … …  
    12961334        }
    12971335}
    1298        
     1336/*}}}*/
     1337/*FUNCTION Penta Misfit {{{1*/
    12991338#undef __FUNCT__
    13001339#define __FUNCT__ "Penta::Misfit"
    13011340double Penta::Misfit(void* inputs,int analysis_type,int sub_analysis_type){
    1302        
     1341
    13031342        double J;
    13041343        Tria* tria=NULL;
    1305        
     1344
    13061345        /*If on water, return 0: */
    13071346        if(onwater)return 0;
    … …  
    13301369        }
    13311370}
    1332                
     1371/*}}}*/
     1372/*FUNCTION Penta SpawnTria {{{1*/
    13331373#undef __FUNCT__
    13341374#define __FUNCT__ "Penta::SpawnTria"
    … …  
    13451385        double   tria_accumulation[3];
    13461386        double   tria_geothermalflux[3];
    1347        
     1387
    13481388        /*configuration: */
    13491389        int tria_node_ids[3];
    … …  
    13701410        tria_melting[1]=melting[g1];
    13711411        tria_melting[2]=melting[g2];
    1372        
     1412
    13731413        tria_accumulation[0]=accumulation[g0];
    13741414        tria_accumulation[1]=accumulation[g1];
    … …  
    14011441
    14021442}
    1403 
    1404 
     1443/*}}}*/
     1444/*FUNCTION Penta GetDofList {{{1*/
    14051445void  Penta::GetDofList(int* doflist,int* pnumberofdofspernode){
    14061446
    … …  
    14081448        int doflist_per_node[MAXDOFSPERNODE];
    14091449        int numberofdofspernode;
    1410        
     1450
    14111451        for(i=0;i<6;i++){
    14121452                nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
    … …  
    14201460
    14211461}
    1422 
     1462/*}}}*/
     1463/*FUNCTION Penta GetDofList1 {{{1*/
    14231464void  Penta::GetDofList1(int* doflist){
    14241465
    … …  
    14291470
    14301471}
     1472/*}}}*/
     1473/*FUNCTION Penta GetStrainRate {{{1*/
    14311474#undef __FUNCT__
    14321475#define __FUNCT__ "Penta::GetStrainRate"
    … …  
    14421485        GetB(&B[0][0], xyz_list, gauss_l1l2l3l4);
    14431486
    1444         #ifdef _ISSM_DEBUG_
     1487#ifdef _ISSM_DEBUG_
    14451488        printf("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
    14461489        printf("              [ %lf   %lf  \n",B[1][0],B[1][1]);
    … …  
    14481491        printf("              [ %lf   %lf  ]\n",B[3][0],B[3][1]);
    14491492        printf("              [ %lf   %lf  ]\n",B[4][0],B[4][1]);
    1450        
     1493
    14511494        printf("B for grid2 : [ %lf   %lf  \n",B[0][2],B[0][3]);
    14521495        printf("              [ %lf   %lf  \n",B[1][2],B[1][3]);
    … …  
    14541497        printf("              [ %lf   %lf  ]\n",B[3][2],B[3][3]);
    14551498        printf("              [ %lf   %lf  ]\n",B[4][2],B[4][3]);
    1456        
     1499
    14571500        printf("B for grid3 : [ %lf   %lf  \n", B[0][4],B[0][5]);
    14581501        printf("              [ %lf   %lf  \n", B[1][4],B[1][5]);
    … …  
    14601503        printf("              [ %lf   %lf  ]\n",B[3][4],B[3][5]);
    14611504        printf("              [ %lf   %lf  ]\n",B[4][4],B[4][5]);
    1462        
     1505
    14631506        printf("B for grid4 : [ %lf   %lf  \n", B[0][6],B[0][7]);
    14641507        printf("              [ %lf   %lf  \n", B[1][6],B[1][7]);
    … …  
    14661509        printf("              [ %lf   %lf  ]\n",B[3][6],B[3][7]);
    14671510        printf("              [ %lf   %lf  ]\n",B[4][6],B[4][7]);
    1468                                
     1511
    14691512        printf("B for grid5 : [ %lf   %lf  \n", B[0][8],B[0][9]);
    14701513        printf("              [ %lf   %lf  \n", B[1][8],B[1][9]);
    … …  
    14791522        printf("              [ %lf   %lf  ]\n",B[4][10],B[4][11]);
    14801523
    1481         #endif
     1524#endif
    14821525
    14831526        /*Multiply B by velocity, to get strain rate: */
    14841527        MatrixMultiply( &B[0][0],5,NDOF2*numgrids,0,
    1485                                       velocity,NDOF2*numgrids,1,0,
    1486                                                   epsilon,0);
    1487 
    1488 }
    1489 
     1528                                velocity,NDOF2*numgrids,1,0,
     1529                                epsilon,0);
     1530
     1531}
     1532/*}}}*/
     1533/*FUNCTION Penta GetB {{{1*/
    14901534#undef __FUNCT__
    14911535#define __FUNCT__ "Penta::GetB"
    … …  
    15041548         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    15051549         */
    1506        
     1550
    15071551        int i;
    15081552        const int numgrids=6;
    … …  
    15151559        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    15161560
    1517         #ifdef _ISSM_DEBUG_
     1561#ifdef _ISSM_DEBUG_
    15181562        for (i=0;i<numgrids;i++){
    15191563                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]);
    15201564        }
    1521         #endif
     1565#endif
    15221566
    15231567        /*Build B: */
    … …  
    15341578                *(B+NDOF2*numgrids*3+NDOF2*i)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
    15351579                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    1536                
     1580
    15371581                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    15381582                *(B+NDOF2*numgrids*4+NDOF2*i+1)=(float).5*dh1dh2dh3dh4dh5dh6_basic[2][i];
    … …  
    15401584
    15411585}
    1542 
    1543 
     1586/*}}}*/
     1587/*FUNCTION Penta GetBPrime {{{1*/
    15441588#undef __FUNCT__
    15451589#define __FUNCT__ "Penta::GetBPrime"
    … …  
    15581602         * We assume B has been allocated already, of size: 5x(NDOF2*numgrids)
    15591603         */
    1560        
     1604
    15611605        int i;
    15621606        const int NDOF3=3;
    … …  
    15691613        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    15701614
    1571         #ifdef _ISSM_DEBUG_
     1615#ifdef _ISSM_DEBUG_
    15721616        for (i=0;i<numgrids;i++){
    15731617                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]);
    15741618        }
    1575         #endif
     1619#endif
    15761620
    15771621        /*Build BPrime: */
    … …  
    15881632                *(B+NDOF2*numgrids*3+NDOF2*i)=dh1dh2dh3dh4dh5dh6_basic[2][i];
    15891633                *(B+NDOF2*numgrids*3+NDOF2*i+1)=0.0;
    1590                
     1634
    15911635                *(B+NDOF2*numgrids*4+NDOF2*i)=0.0;
    15921636                *(B+NDOF2*numgrids*4+NDOF2*i+1)=dh1dh2dh3dh4dh5dh6_basic[2][i];
    15931637        }
    15941638}
    1595 
     1639/*}}}*/
     1640/*FUNCTION Penta GetJacobianDeterminant {{{1*/
    15961641#undef __FUNCT__
    15971642#define __FUNCT__ "Penta::GetJacobianDeterminant"
    … …  
    16011646         * the determinant of it: */
    16021647        const int NDOF3=3;
    1603        
     1648
    16041649        double J[NDOF3][NDOF3];
    16051650
    … …  
    16111656        }
    16121657}
    1613 
     1658/*}}}*/
     1659/*FUNCTION Penta GetNodalFunctionsDerivativesBasic {{{1*/
    16141660#undef __FUNCT__
    16151661#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasic"
    16161662void Penta::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3dh4dh5dh6_basic,double* xyz_list, double* gauss_l1l2l3l4){
    1617        
     1663
    16181664        /*This routine returns the values of the nodal functions derivatives  (with respect to the basic coordinate system: */
    16191665
    1620        
     1666
    16211667        int i;
    16221668        const int NDOF3=3;
    … …  
    16461692
    16471693}
    1648 
     1694/*}}}*/
     1695/*FUNCTION Penta GetJacobian {{{1*/
    16491696#undef __FUNCT__
    16501697#define __FUNCT__ "Penta::GetJacobian"
    … …  
    16531700        const int NDOF3=3;
    16541701        int i,j;
    1655        
     1702
    16561703        /*The Jacobian is constant over the element, discard the gaussian points.
    16571704         * J is assumed to have been allocated of size NDOF2xNDOF2.*/
    … …  
    16631710        double y1,y2,y3,y4,y5,y6;
    16641711        double z1,z2,z3,z4,z5,z6;
    1665        
     1712
    16661713        double sqrt3=sqrt(3.0);
    1667        
     1714
    16681715        /*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
    16691716        A1=gauss_l1l2l3l4[0];
    … …  
    16811728        x5=*(xyz_list+3*4+0);
    16821729        x6=*(xyz_list+3*5+0);
    1683        
     1730
    16841731        y1=*(xyz_list+3*0+1);
    16851732        y2=*(xyz_list+3*1+1);
    … …  
    17091756        *(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);
    17101757
    1711         #ifdef _ISSM_DEBUG_
     1758#ifdef _ISSM_DEBUG_
    17121759        for(i=0;i<3;i++){
    17131760                for (j=0;j<3;j++){
    … …  
    17161763                printf("\n");
    17171764        }
    1718         #endif
    1719 }
    1720 
     1765#endif
     1766}
     1767/*}}}*/
     1768/*FUNCTION Penta GetNodalFunctionsDerivativesParams {{{1*/
    17211769#undef __FUNCT__
    17221770#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParams"
    17231771void Penta::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3dl4dl5dl6,double* gauss_l1l2l3l4){
    1724        
     1772
    17251773        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    17261774         * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */
    … …  
    17291777        double A1,A2,A3,z;
    17301778        double sqrt3=sqrt(3.0);
    1731        
     1779
    17321780        A1=gauss_l1l2l3l4[0]; //first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*sqrt(3));
    17331781        A2=gauss_l1l2l3l4[1]; //second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*sqrt(3));
    … …  
    17661814        *(dl1dl2dl3dl4dl5dl6+numgrids*2+5)=1.0/2.0*A3;
    17671815}
    1768 
     1816/*}}}*/
     1817/*FUNCTION Penta GetJacobianInvert {{{1*/
    17691818#undef __FUNCT__
    17701819#define __FUNCT__ "Penta::GetJacobianInvert"
    … …  
    17801829        MatrixInverse(Jinv,NDOF3,NDOF3,NULL,0,&Jdet);
    17811830}
    1782 
    1783 
    1784 
     1831/*}}}*/
     1832/*FUNCTION Penta CreatePVectorDiagnosticHoriz {{{1*/
    17851833#undef __FUNCT__
    17861834#define __FUNCT__ "Penta::CreatePVectorDiagnosticHoriz"
    … …  
    17961844        int          doflist[numdof];
    17971845        int          numberofdofspernode;
    1798        
     1846
    17991847        /* parameters: */
    18001848        double  slope[NDOF2];
    … …  
    18741922
    18751923                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);
    1876                 #ifdef _ISSM_DEBUG_
     1924#ifdef _ISSM_DEBUG_
    18771925                for (i=0;i<num_area_gauss;i++){
    18781926                        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));
    … …  
    18811929                        printf("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
    18821930                }
    1883                 #endif
    1884        
     1931#endif
     1932
    18851933                /* Start  looping on the number of gaussian points: */
    18861934                for (ig1=0; ig1<num_area_gauss; ig1++){
    18871935                        for (ig2=0; ig2<num_vert_gauss; ig2++){
    1888                        
     1936
    18891937                                /*Pick up the gaussian point: */
    18901938                                gauss_weight1=*(area_gauss_weights+ig1);
    18911939                                gauss_weight2=*(vert_gauss_weights+ig2);
    18921940                                gauss_weight=gauss_weight1*gauss_weight2;
    1893                                
     1941
    18941942                                gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    18951943                                gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    18961944                                gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    18971945                                gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
    1898                
     1946
    18991947                                /*Compute thickness at gaussian point: */
    19001948                                GetParameterValue(&thickness, &h[0],gauss_l1l2l3l4);
    … …  
    19051953                                /* Get Jacobian determinant: */
    19061954                                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
    1907                
    1908                                  /*Get nodal functions: */
     1955
     1956                                /*Get nodal functions: */
    19091957                                GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
    19101958
    … …  
    19301978        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    19311979
    1932         cleanup_and_return:
     1980cleanup_and_return:
    19331981        xfree((void**)&first_gauss_area_coord);
    19341982        xfree((void**)&second_gauss_area_coord);
    … …  
    19391987
    19401988}
    1941 
    1942 
     1989/*}}}*/
     1990/*FUNCTION Penta GetParameterValue {{{1*/
    19431991#undef __FUNCT__
    19441992#define __FUNCT__ "Penta::GetParameterValue"
    19451993void Penta::GetParameterValue(double* pvalue, double* v_list,double* gauss_l1l2l3l4){
    1946        
     1994
    19471995        const int numgrids=6;
    19481996        double l1l2l3l4l5l6[numgrids];
    … …  
    19522000        *pvalue=l1l2l3l4l5l6[0]*v_list[0]+l1l2l3l4l5l6[1]*v_list[1]+l1l2l3l4l5l6[2]*v_list[2]+l1l2l3l4l5l6[3]*v_list[3]+l1l2l3l4l5l6[4]*v_list[4]+l1l2l3l4l5l6[5]*v_list[5];
    19532001}
    1954 
     2002/*}}}*/
     2003/*FUNCTION Penta GetParameterDerivativeValue {{{1*/
    19552004#undef __FUNCT__
    19562005#define __FUNCT__ "Penta::GetParameterDerivativeValue"
    19572006void Penta::GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_l1l2l3l4){
    1958                                
     2007
    19592008        /*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:
    19602009         *   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;
    … …  
    19712020        /*Get dh1dh2dh3dh4dh5dh6_basic in basic coordinate system: */
    19722021        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    1973        
     2022
    19742023        *(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];
    1975 ;
     2024        ;
    19762025        *(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];
    19772026
    … …  
    19792028
    19802029}
    1981 
     2030/*}}}*/
     2031/*FUNCTION Penta GetNodalFunctions {{{1*/
    19822032#undef __FUNCT__
    19832033#define __FUNCT__ "Penta::GetNodalFunctions"
    19842034void Penta::GetNodalFunctions(double* l1l2l3l4l5l6, double* gauss_l1l2l3l4){
    1985        
     2035
    19862036        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    19872037
    … …  
    19892039
    19902040        l1l2l3l4l5l6[1]=gauss_l1l2l3l4[1]*(1-gauss_l1l2l3l4[3])/2.0;
    1991        
     2041
    19922042        l1l2l3l4l5l6[2]=gauss_l1l2l3l4[2]*(1-gauss_l1l2l3l4[3])/2.0;
    19932043
    … …  
    19952045
    19962046        l1l2l3l4l5l6[4]=gauss_l1l2l3l4[1]*(1+gauss_l1l2l3l4[3])/2.0;
    1997        
     2047
    19982048        l1l2l3l4l5l6[5]=gauss_l1l2l3l4[2]*(1+gauss_l1l2l3l4[3])/2.0;
    19992049
    20002050}
    2001 
     2051/*}}}*/
     2052/*FUNCTION Penta FieldExtrude {{{1*/
    20022053#undef __FUNCT__
    20032054#define __FUNCT__ "Penta::FieldExtrude"
    20042055void  Penta::FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed){
    2005                
     2056
    20062057        /* node data: */
    20072058        const int    numgrids=6;
    … …  
    20442095                         * inserting the same field value into field, until we reach the surface: */
    20452096                        for(i=0;i<3;i++){
    2046                
     2097
    20472098                                node=nodes[i]; //base nodes
    2048                
     2099
    20492100                                /*get field for this base node: */
    20502101                                fieldel[0]=field_serial[doflist[numberofdofspernode*i+0]];
    … …  
    20962147                }
    20972148                else if (
    2098                                 (strcmp(field_name,"thickness")==0) ||
    2099                                 (strcmp(field_name,"surface")==0)  ||
    2100                                 (strcmp(field_name,"bed")==0)  ||
    2101                                 (strcmp(field_name,"slopex")==0)  ||
    2102                                 (strcmp(field_name,"slopey")==0)
    2103                                 ){
     2149                                        (strcmp(field_name,"thickness")==0) ||
     2150                                        (strcmp(field_name,"surface")==0)  ||
     2151                                        (strcmp(field_name,"bed")==0)  ||
     2152                                        (strcmp(field_name,"slopex")==0)  ||
     2153                                        (strcmp(field_name,"slopey")==0)
     2154                                  ){
    21042155
    21052156                        /* node data: */
    … …  
    21082159                        int          nodedofs;
    21092160                        double       fieldel;
    2110                                
     2161
    21112162                        GetDofList(&doflist[0],&numberofdofspernode);
    21122163
    … …  
    21152166                         * inserting the same thickness value into tg, until we reach the surface: */
    21162167                        for(i=0;i<3;i++){
    2117                
     2168
    21182169                                node=nodes[i]; //base nodes
    2119                
     2170
    21202171                                /*get velocity for this base node: */
    21212172                                fieldel=field_serial[doflist[numberofdofspernode*i+0]];
    … …  
    21392190        } //if (extrude)
    21402191}
    2141 
     2192/*}}}*/
     2193/*FUNCTION Penta GetB_vert {{{1*/
    21422194#undef __FUNCT__
    21432195#define __FUNCT__ "Penta:GetB_vert"
    … …  
    21462198
    21472199        /*      Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
    2148         where hi is the interpolation function for grid i.*/
     2200                where hi is the interpolation function for grid i.*/
    21492201
    21502202        int i;
    … …  
    21562208        GetNodalFunctionsDerivativesBasic(&dh1dh2dh3dh4dh5dh6_basic[0][0],xyz_list, gauss_l1l2l3l4);
    21572209
    2158         #ifdef _ISSM_DEBUG_
     2210#ifdef _ISSM_DEBUG_
    21592211        for (i=0;i<numgrids;i++){
    21602212                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]);
    21612213        }
    2162         #endif
     2214#endif
    21632215
    21642216        /*Build B: */
    … …  
    21662218                B[i]=dh1dh2dh3dh4dh5dh6_basic[2][i]; 
    21672219        }
    2168        
    2169 }
    2170 
    2171 
     2220
     2221}
     2222/*}}}*/
     2223/*FUNCTION Penta GetBPrime_vert {{{1*/
    21722224#undef __FUNCT__
    21732225#define __FUNCT__ "Penta:GetBPrime_vert"
    … …  
    21772229
    21782230        int i;
    2179                                
     2231
    21802232        GetNodalFunctions(B, gauss_l1l2l3l4);
    21812233
    21822234}
    2183 
     2235/*}}}*/
     2236/*FUNCTION Penta CreatePVectorDiagnosticVert {{{1*/
    21842237#undef __FUNCT__
    21852238#define __FUNCT__ "Penta:CreatePVectorDiagnosticVert"
    … …  
    21962249        int          doflist[numdof];
    21972250        int          numberofdofspernode;
    2198        
     2251
    21992252        /* gaussian points: */
    22002253        int     num_gauss,ig;
    … …  
    22592312        /* recover input parameters: */
    22602313        if(!inputs->Recover("velocity",&vx_list[0],1,dofs1,numgrids,(void**)nodes))throw ErrorException(__FUNCT__," cannot compute vertical velocity without horizontal velocity!");
    2261             inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
     2314        inputs->Recover("velocity",&vy_list[0],1,dofs2,numgrids,(void**)nodes);
    22622315
    22632316        /*Get gaussian points and weights :*/
    … …  
    22662319
    22672320        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);
    2268         #ifdef _ISSM_DEBUG_
     2321#ifdef _ISSM_DEBUG_
    22692322        for (i=0;i<num_area_gauss;i++){
    22702323                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));
    … …  
    22732326                printf("Vert Gauss coord %i: %lf Weight: %lf\n",i,*(fourth_gauss_vert_coord+i),*(vert_gauss_weights+i));
    22742327        }
    2275         #endif
     2328#endif
    22762329
    22772330        /* Start  looping on the number of gaussian points: */
    22782331        for (ig1=0; ig1<num_area_gauss; ig1++){
    22792332                for (ig2=0; ig2<num_vert_gauss; ig2++){
    2280                
     2333
    22812334                        /*Pick up the gaussian point: */
    22822335                        gauss_weight1=*(area_gauss_weights+ig1);
    22832336                        gauss_weight2=*(vert_gauss_weights+ig2);
    22842337                        gauss_weight=gauss_weight1*gauss_weight2;
    2285                        
     2338
    22862339                        gauss_l1l2l3l4[0]=*(first_gauss_area_coord+ig1);
    22872340                        gauss_l1l2l3l4[1]=*(second_gauss_area_coord+ig1);
    22882341                        gauss_l1l2l3l4[2]=*(third_gauss_area_coord+ig1);
    22892342                        gauss_l1l2l3l4[3]=*(fourth_gauss_vert_coord+ig2);
    2290        
     2343
    22912344                        /*Get velocity derivative, with respect to x and y: */
    22922345                        GetParameterDerivativeValue(&du[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3l4);
    … …  
    22942347                        dudx=du[0];
    22952348                        dvdy=dv[1];
    2296                        
     2349
    22972350
    22982351                        /* Get Jacobian determinant: */
    22992352                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3l4);
    2300                         #ifdef _ISSM_DEBUG_
     2353#ifdef _ISSM_DEBUG_
    23012354                        printf("Element id %i Jacobian determinant: %lf\n",GetId(),Jdet);
    2302                         #endif
    2303                
    2304                          /*Get nodal functions: */
     2355#endif
     2356
     2357                        /*Get nodal functions: */
    23052358                        GetNodalFunctions(l1l2l3l4l5l6, gauss_l1l2l3l4);
    23062359
    … …  
    23192372        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    23202373
    2321         cleanup_and_return:
     2374cleanup_and_return:
    23222375        xfree((void**)&first_gauss_area_coord);
    23232376        xfree((void**)&second_gauss_area_coord);
    … …  
    23262379        xfree((void**)&area_gauss_weights);
    23272380        xfree((void**)&vert_gauss_weights);
    2328 
    2329 
    2330 }
    2331 
     2381}
     2382/*}}}*/
     2383/*FUNCTION Penta ComputePressure {{{1*/
    23322384#undef __FUNCT__
    23332385#define __FUNCT__ "Penta::ComputePressure"
    … …  
    23432395        /*If on water, skip: */
    23442396        if(onwater)return;
    2345                
     2397
    23462398        /*Get node data: */
    23472399        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    2348        
     2400
    23492401        /*pressure is lithostatic: */
    23502402        //md.pressure=md.rho_ice*md.g*(md.surface-md.z); a la matlab
    … …  
    23592411                pressure[i]=rho_ice*g*(s[i]-xyz_list[i][2]);
    23602412        }
    2361        
     2413
    23622414        /*plug local pressure values into global pressure vector: */
    23632415        VecSetValues(pg,numgrids,doflist,(const double*)pressure,INSERT_VALUES);
    23642416
    23652417}
    2366 
    2367 
     2418/*}}}*/
     2419/*FUNCTION Penta CreateKMatrixSlopeCompute {{{1*/
    23682420#undef __FUNCT__
    23692421#define __FUNCT__ "Penta::CreateKMatrixSlopeCompute"
    … …  
    23762428        /*If on water, skip: */
    23772429        if(onwater)return;
    2378        
     2430
    23792431        /*Is this element on the bed? :*/
    23802432        if(!onbed)return;
    … …  
    23872439
    23882440}
    2389 
     2441/*}}}*/
     2442/*FUNCTION Penta CreatePVectorSlopeCompute {{{1*/
    23902443#undef __FUNCT__
    23912444#define __FUNCT__ "Penta::CreatePVectorSlopeCompute"
    23922445
    23932446void Penta::CreatePVectorSlopeCompute( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    2394        
     2447
    23952448        /*Collapsed formulation: */
    23962449        Tria*  tria=NULL;
    … …  
    23982451        /*If on water, skip: */
    23992452        if(onwater)return;
    2400        
     2453
    24012454        /*Is this element on the bed? :*/
    24022455        if(!onbed)return;
    … …  
    24082461        return;
    24092462}
    2410 
     2463/*}}}*/
     2464/*FUNCTION Penta CreateKMatrixPrognostic {{{1*/
    24112465#undef __FUNCT__
    24122466#define __FUNCT__ "Penta::CreateKMatrixPrognostic"
    … …  
    24192473        /*If on water, skip: */
    24202474        if(onwater)return;
    2421        
     2475
    24222476        /*Is this element on the bed? :*/
    24232477        if(!onbed)return;
    … …  
    24302484
    24312485}
    2432 
     2486/*}}}*/
     2487/*FUNCTION Penta CreatePVectorPrognostic {{{1*/
    24332488#undef __FUNCT__
    24342489#define __FUNCT__ "Penta::CreatePVectorPrognostic"
    24352490
    24362491void Penta::CreatePVectorPrognostic( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){
    2437        
     2492
    24382493        /*Collapsed formulation: */
    24392494        Tria*  tria=NULL;
    … …  
    24412496        /*If on water, skip: */
    24422497        if(onwater)return;
    2443        
     2498
    24442499        /*Is this element on the bed? :*/
    24452500        if(!onbed)return;
    … …  
    24512506        return;
    24522507}
    2453 
     2508/*}}}*/
     2509/*FUNCTION Penta ReduceMatrixStokes {{{1*/
    24542510#undef __FUNCT__
    24552511#define __FUNCT__ "ReduceMatrixStokes"
    24562512void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
    2457                
     2513
    24582514        int i,j;
    24592515
    … …  
    24972553                }
    24982554        }
    2499 
    2500 
    2501 }
    2502 
     2555}
     2556/*}}}*/
     2557/*FUNCTION Penta GetMatrixInvert {{{1*/
    25032558#undef __FUNCT__
    25042559#define __FUNCT__ "GetMatrixInvert"
    … …  
    25092564        double det;
    25102565        int calculationdof=3;
    2511        
     2566
    25122567        /*Take the matrix components: */
    25132568        a=*(Ke+calculationdof*0+0);
    … …  
    25222577
    25232578        det=a*(e*i-f*h)-b*(d*i-f*g)+c*(d*h-e*g);
    2524        
     2579
    25252580        *(Ke_invert+calculationdof*0+0)=(e*i-f*h)/det;
    25262581        *(Ke_invert+calculationdof*0+1)=(c*h-b*i)/det;
    … …  
    25342589
    25352590}
    2536 
     2591/*}}}*/
     2592/*FUNCTION Penta SurfaceNormal {{{1*/
    25372593#undef __FUNCT__
    25382594#define __FUNCT__ "Penta::SurfaceNormal"
    2539 
    25402595void Penta::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){
    25412596
    … …  
    25622617
    25632618}
    2564 
     2619/*}}}*/
     2620/*FUNCTION Penta GetStrainRateStokes {{{1*/
    25652621#undef __FUNCT__
    25662622#define __FUNCT__ "Penta::GetStrainRateStokes"
    … …  
    26012657        MatrixMultiply( &B_reduced[0][0],6,DOFVELOCITY*numgrids, 0, velocity,DOFVELOCITY*numgrids,1,0,epsilon, 0);
    26022658}
    2603 
     2659/*}}}*/
     2660/*FUNCTION Penta GetBStokes {{{1*/
    26042661#undef __FUNCT__
    26052662#define __FUNCT__ "Penta::GetBStokes"
    … …  
    26092666         * For grid i, Bi can be expressed in the basic coordinate system
    26102667         * by:                          Bi_basic=[ dh/dx          0             0       0  ]
    2611                 *                                       [   0           dh/dy           0       0  ]
    2612                 *                                       [   0             0           dh/dy     0  ]
    2613                 *                                       [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
    2614                 *                                       [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
    2615                 *                                       [   0          1/2*dh/dz    1/2*dh/dy   0  ]
    2616                 *                                       [   0             0             0       h  ]
    2617                 *                                       [ dh/dx         dh/dy         dh/dz     0  ]
    2618                 *       where h is the interpolation function for grid i.
    2619                 *       Same thing for Bb except the last column that does not exist.
    2620                 */
    2621        
     2668         *                                      [   0           dh/dy           0       0  ]
     2669         *                                      [   0             0           dh/dy     0  ]
     2670         *                                      [ 1/2*dh/dy    1/2*dh/dx        0       0  ]
     2671         *                                      [ 1/2*dh/dz       0         1/2*dh/dx   0  ]
     2672         *                                      [   0          1/2*dh/dz    1/2*dh/dy   0  ]
     2673         *                                      [   0             0             0       h  ]
     2674         *                                      [ dh/dx         dh/dy         dh/dz     0  ]
     2675         *      where h is the interpolation function for grid i.
     2676         *      Same thing for Bb except the last column that does not exist.
     2677         */
     2678
    26222679        int i;
    26232680        const int calculationdof=3;
    … …  
    26332690
    26342691        GetNodalFunctions(l1l6, gauss_coord);
    2635        
    2636         #ifdef _ISSM_DEBUG_
     2692
     2693#ifdef _ISSM_DEBUG_
    26372694        for (i=0;i<7;i++){
    26382695                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]);
    26392696        }
    26402697
    2641         #endif
     2698#endif
    26422699
    26432700        /*Build B: */
    … …  
    26812738
    26822739}
    2683 
     2740/*}}}*/
     2741/*FUNCTION Penta GetBprimeStokes {{{1*/
    26842742#undef __FUNCT__
    26852743#define __FUNCT__ "Penta::GetBprimeStokes"
    … …  
    26872745
    26882746        /*      Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2.
    2689         *       For grid i, Bi' can be expressed in the basic coordinate system
    2690         *       by:
    2691         *                               Bi_basic'=[ dh/dx   0          0       0]
    2692         *                                        [   0      dh/dy      0       0]
    2693         *                                        [   0      0         dh/dz    0]
    2694         *                                        [  dh/dy   dh/dx      0       0]
    2695         *                                        [  dh/dz   0        dh/dx     0]
    2696         *                                        [   0      dh/dz    dh/dy     0]
    2697         *                                        [  dh/dx   dh/dy    dh/dz     0]
    2698         *                                        [   0      0          0       h]
    2699         *       where h is the interpolation function for grid i.
    2700         *
    2701         *       Same thing for the bubble fonction except that there is no fourth column
    2702         */
     2747         *      For grid i, Bi' can be expressed in the basic coordinate system
     2748         *      by:
     2749         *                              Bi_basic'=[ dh/dx   0          0       0]
     2750         *                                       [   0      dh/dy      0       0]
     2751         *                                       [   0      0         dh/dz    0]
     2752         *                                       [  dh/dy   dh/dx      0       0]
     2753         *                                       [  dh/dz   0        dh/dx     0]
     2754         *                                       [   0      dh/dz    dh/dy     0]
     2755         *                                       [  dh/dx   dh/dy    dh/dz     0]
     2756         *                                       [   0      0          0       h]
     2757         *      where h is the interpolation function for grid i.
     2758         *
     2759         *      Same thing for the bubble fonction except that there is no fourth column
     2760         */
    27032761
    27042762        int i;
    … …  
    27102768        double l1l6[numgrids];
    27112769
    2712 
    27132770        /*Get dh1dh7 in basic coordinate system: */
    27142771        GetNodalFunctionsDerivativesBasicStokes(&dh1dh7_basic[0][0],xyz_list, gauss_coord);
    27152772
    27162773        GetNodalFunctions(l1l6, gauss_coord);
    2717        
    2718         #ifdef _ISSM_DEBUG_
     2774
     2775#ifdef _ISSM_DEBUG_
    27192776        for (i=0;i<6;i++){
    27202777                printf("Node %i  dh/dx=%lf dh/dy=%lf dh/dz=%lf \n",i,dh1dh7_basic[0][i],dh1dh7_basic[1][i]);
    27212778        }
    27222779
    2723         #endif
     2780#endif
    27242781
    27252782        /*B_primeuild B_prime: */
    … …  
    27632820
    27642821}
     2822/*}}}*/
     2823/*FUNCTION Penta GetLStokes {{{1*/
    27652824#undef __FUNCT__
    27662825#define __FUNCT__ "Penta::GetLStokes"
    … …  
    27682827
    27692828        /*
    2770         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
    2771         * For grid i, Li can be expressed in the basic coordinate system
    2772         * by:
    2773         *       Li_basic=[ h    0    0   0]
    2774         *                [ 0    h    0   0]
    2775         *                [ 0    0    h   0]
    2776         *                [ 0    0    h   0]
    2777         *                [ h    0    0   0]
    2778         *                [ 0    h    0   0]
    2779         *                [ h    0    0   0]
    2780         *                [ 0    h    0   0]
    2781         *                [ 0    0    h   0]
    2782         *                [ 0    0    h   0]
    2783         *                [ 0    0    h   0]
    2784         *                [ h    0    0   0]
    2785         *                [ 0    h    0   0]
    2786         *                [ 0    0    h   0]
    2787         * where h is the interpolation function for grid i.
    2788         */
     2829         * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof.
     2830         * For grid i, Li can be expressed in the basic coordinate system
     2831         * by:
     2832         *       Li_basic=[ h    0    0   0]
     2833         *               [ 0    h    0   0]
     2834         *               [ 0    0    h   0]
     2835         *               [ 0    0    h   0]
     2836         *               [ h    0    0   0]
     2837         *               [ 0    h    0   0]
     2838         *               [ h    0    0   0]
     2839         *               [ 0    h    0   0]
     2840         *               [ 0    0    h   0]
     2841         *               [ 0    0    h   0]
     2842         *               [ 0    0    h   0]
     2843         *               [ h    0    0   0]
     2844         *               [ 0    h    0   0]
     2845         *               [ 0    0    h   0]
     2846         * where h is the interpolation function for grid i.
     2847         */
    27892848
    27902849        int i;
    … …  
    27992858        l1l2l3[1]=gauss_coord_tria[1];
    28002859        l1l2l3[2]=gauss_coord_tria[2];
    2801        
    2802 
    2803         #ifdef _DELUG_
     2860
     2861
     2862#ifdef _DELUG_
    28042863        for (i=0;i<3;i++){
    28052864                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
    28062865        }
    2807         #endif
     2866#endif
    28082867
    28092868        /*Build LStokes: */
    … …  
    28652924                *(LStokes+num_dof*numgrids2d*13+num_dof*i+2)=l1l2l3[i];
    28662925                *(LStokes+num_dof*numgrids2d*13+num_dof*i+3)=0;
    2867        
    2868         }
    2869 }
    2870 
    2871 
     2926
     2927        }
     2928}
     2929/*}}}*/
     2930/*FUNCTION Penta GetLprimeStokes {{{1*/
    28722931#undef __FUNCT__
    28732932#define __FUNCT__ "Penta::GetLprimeStokes"
    … …  
    28752934
    28762935        /*
    2877         * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
    2878         * For grid i, Lpi can be expressed in the basic coordinate system
    2879         * by:
    2880         *       Lpi_basic=[ h    0    0   0]
    2881         *                [ 0    h    0   0]
    2882         *                [ h    0    0   0]
    2883         *                [ 0    h    0   0]
    2884         *                [ 0    0    h   0]
    2885         *                [ 0    0    h   0]
    2886         *                [ 0    0  dh/dz 0]
    2887         *                [ 0    0  dh/dz 0]
    2888         *                [ 0    0  dh/dz 0]
    2889         *                [dh/dz 0  dh/dx 0]
    2890         *                [ 0 dh/dz dh/dy 0]
    2891         *                [ 0    0    0   h]
    2892         *                [ 0    0    0   h]
    2893         *                [ 0    0    0   h]
    2894         * where h is the interpolation function for grid i.
    2895         */
     2936         * Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof.
     2937         * For grid i, Lpi can be expressed in the basic coordinate system
     2938         * by:
     2939         *       Lpi_basic=[ h    0    0   0]
     2940         *               [ 0    h    0   0]
     2941         *               [ h    0    0   0]
     2942         *               [ 0    h    0   0]
     2943         *               [ 0    0    h   0]
     2944         *               [ 0    0    h   0]
     2945         *               [ 0    0  dh/dz 0]
     2946         *               [ 0    0  dh/dz 0]
     2947         *               [ 0    0  dh/dz 0]
     2948         *               [dh/dz 0  dh/dx 0]
     2949         *               [ 0 dh/dz dh/dy 0]
     2950         *               [ 0    0    0   h]
     2951         *               [ 0    0    0   h]
     2952         *               [ 0    0    0   h]
     2953         * where h is the interpolation function for grid i.
     2954         */
    28962955
    28972956        int i;
    … …  
    29072966        l1l2l3[1]=gauss_coord_tria[1];
    29082967        l1l2l3[2]=gauss_coord_tria[2];
    2909        
     2968
    29102969        GetNodalFunctionsDerivativesBasic(&dh1dh6_basic[0][0],xyz_list,gauss_coord);
    29112970
    2912         #ifdef _DELUG_
     2971#ifdef _DELUG_
    29132972        for (i=0;i<3;i++){
    29142973                printf("Node %i  h=%lf \n",i,l1l2l3[i]);
    29152974        }
    2916         #endif
     2975#endif
    29172976
    29182977        /*Build LprimeStokes: */
    … …  
    29743033                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+2)=0;
    29753034                *(LprimeStokes+num_dof*numgrids2d*13+num_dof*i+3)=l1l2l3[i];
    2976        
    2977         }
    2978        
    2979 }
    2980 
     3035
     3036        }
     3037}
     3038/*}}}*/
     3039/*FUNCTION Penta GetNodalFunctionsDerivativesBasicStokes {{{1*/
    29813040#undef __FUNCT__
    29823041#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesBasicStokes"
    29833042void Penta::GetNodalFunctionsDerivativesBasicStokes(double* dh1dh7_basic,double* xyz_list, double* gauss_coord){
    2984        
     3043
    29853044        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    29863045         * basic coordinate system: */
    … …  
    30133072
    30143073}
    3015 
    3016 
     3074/*}}}*/
     3075/*FUNCTION Penta GetNodalFunctionsDerivativesParamsStokes {{{1*/
    30173076#undef __FUNCT__
    30183077#define __FUNCT__ "Penta::GetNodalFunctionsDerivativesParamsStokes"
    30193078void Penta::GetNodalFunctionsDerivativesParamsStokes(double* dl1dl7,double* gauss_coord){
    3020        
     3079
    30213080        /*This routine returns the values of the nodal functions derivatives  (with respect to the
    30223081         * natural coordinate system) at the gaussian point. */
    … …  
    30653124
    30663125}
    3067 
     3126/*}}}*/
     3127/*FUNCTION Penta CreatePVectorDiagnosticStokes {{{1*/
    30683128#undef __FUNCT__
    30693129#define __FUNCT__ "Penta::CreatePVectorDiagnosticStokes"
    … …  
    31093169        int     area_order=5;
    31103170        int       num_vert_gauss=5;
    3111        
     3171
    31123172        double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    31133173        double  viscosity;
    … …  
    31323192        double     tBD[27][8];
    31333193        double     P_terms[numdof];
    3134        
     3194
    31353195        ParameterInputs* inputs=NULL;
    31363196        Tria*            tria=NULL;
    … …  
    31603220
    31613221        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    3162            get tria gaussian points as well as segment gaussian points. For tria gaussian
    3163            points, order of integration is 2, because we need to integrate the product tB*D*B'
    3164            which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    3165            points, same deal, which yields 3 gaussian points.*/
     3222                get tria gaussian points as well as segment gaussian points. For tria gaussian
     3223                points, order of integration is 2, because we need to integrate the product tB*D*B'
     3224                which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3225                points, same deal, which yields 3 gaussian points.*/
    31663226
    31673227        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    … …  
    31833243                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    31843244                        matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3185                
     3245
    31863246                        /* Get Jacobian determinant: */
    31873247                        GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_coord);
    … …  
    32033263                        GetBStokes(&B[0][0],&xyz_list[0][0],gauss_coord);
    32043264                        GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss_coord);
    3205                        
     3265
    32063266                        /*Get bubble part of Bprime */
    32073267                        for(i=0;i<8;i++){
    … …  
    32423302                        }
    32433303                }
    3244        
     3304
    32453305                xfree((void**)&first_gauss_area_coord); xfree((void**)&second_gauss_area_coord); xfree((void**)&third_gauss_area_coord); xfree((void**)&area_gauss_weights);
    32463306                GaussTria (&num_area_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &area_gauss_weights, 2);
    … …  
    32533313                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
    32543314                        gauss_coord[3]=-1;
    3255                        
     3315
    32563316                        gauss_coord_tria[0]=*(first_gauss_area_coord+igarea);
    32573317                        gauss_coord_tria[1]=*(second_gauss_area_coord+igarea);
    … …  
    32603320                        /*Get the Jacobian determinant */
    32613321                        tria->GetJacobianDeterminant3d(&Jdet2d, &xyz_list_tria[0][0], gauss_coord_tria);
    3262                        
     3322
    32633323                        /* Get bed at gaussian point */
    32643324                        GetParameterValue(&bed,&b[0],gauss_coord);
    … …  
    32663326                        /*Get L matrix : */
    32673327                        tria->GetL(&L[0], &xyz_list[0][0], gauss_coord_tria,1);
    3268                                
     3328
    32693329                        /*Get water_pressure at gaussian point */
    32703330                        water_pressure=gravity*rho_water*bed;
    … …  
    32723332                        /*Get normal vecyor to the bed */
    32733333                        SurfaceNormal(&surface_normal[0],xyz_list_tria);
    3274                        
     3334
    32753335                        bed_normal[0]=-surface_normal[0]; //Program is for surface, so the normal to the bed is the opposite of the result
    32763336                        bed_normal[1]=-surface_normal[1];
    … …  
    32843344                }
    32853345        } //if ( (onbed==1) && (shelf==1))
    3286        
     3346
    32873347        /*Reduce the matrix */
    32883348        ReduceVectorStokes(&Pe_reduced[0], &Ke_temp[0][0], &Pe_temp[0]);
    … …  
    33043364
    33053365}
    3306 
     3366/*}}}*/
     3367/*FUNCTION Penta ReduceVectorStokes {{{1*/
    33073368#undef __FUNCT__
    33083369#define __FUNCT__ "Penta::ReduceVectorStokes"
    33093370void Penta::ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp){
    3310                
     3371
    33113372        int i,j;
    33123373
    … …  
    33473408        }
    33483409}
    3349 
     3410/*}}}*/
     3411/*FUNCTION Penta GetNodalFunctionsStokes {{{1*/
    33503412#undef __FUNCT__
    33513413#define __FUNCT__ "Penta::GetNodalFunctionsStokes"
    33523414void Penta::GetNodalFunctionsStokes(double* l1l7, double* gauss_coord){
    3353        
     3415
    33543416        /*This routine returns the values of the nodal functions  at the gaussian point.*/
    33553417
    … …  
    33763438
    33773439}
    3378 
     3440/*}}}*/
     3441/*FUNCTION Penta CreateKMatrixThermal {{{1*/
    33793442#undef __FUNCT__
    33803443#define __FUNCT__ "Penta::CreateKMatrixThermal"
    33813444void  Penta::CreateKMatrixThermal(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    3382        
     3445
    33833446        /* local declarations */
    33843447        int i,j;
    … …  
    34663529        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    34673530        GetDofList(&doflist[0],&numberofdofspernode);
    3468        
     3531
    34693532        // /*recovre material parameters: */
    34703533        rho_water=matpar->GetRhoWater();
    … …  
    34733536        heatcapacity=matpar->GetHeatCapacity();
    34743537        thermalconductivity=matpar->GetThermalConductivity();
    3475                
     3538
    34763539        /*recover extra inputs from users, dt and velocity: */
    34773540        found=inputs->Recover("velocity",&vxvyvz_list[0][0],3,dofs,numgrids,(void**)nodes);
    … …  
    34883551
    34893552        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    3490            get tria gaussian points as well as segment gaussian points. For tria gaussian
    3491            points, order of integration is 2, because we need to integrate the product tB*D*B'
    3492            which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    3493            points, same deal, which yields 3 gaussian points.: */
    3494        
     3553                get tria gaussian points as well as segment gaussian points. For tria gaussian
     3554                points, order of integration is 2, because we need to integrate the product tB*D*B'
     3555                which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3556                points, same deal, which yields 3 gaussian points.: */
     3557
    34953558        /*Get gaussian points: */
    34963559        area_order=2;
    … …  
    36013664                                }
    36023665                        }
    3603                                        
     3666
    36043667                        /*Add Ke_gaussian to pKe: */
    36053668                        for(i=0;i<numdof;i++){
    … …  
    36153678        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)K_terms,ADD_VALUES);
    36163679
    3617         cleanup_and_return:
     3680cleanup_and_return:
    36183681        xfree((void**)&first_gauss_area_coord);
    36193682        xfree((void**)&second_gauss_area_coord);
    … …  
    36223685        xfree((void**)&vert_gauss_weights);
    36233686        xfree((void**)&vert_gauss_coord);
    3624        
     3687
    36253688        //Ice/ocean heat exchange flux on ice shelf base
    36263689        if(onbed && shelf){
    … …  
    36313694        }
    36323695}
    3633 
     3696/*}}}*/
     3697/*FUNCTION Penta GetB_conduct {{{1*/
    36343698#undef __FUNCT__
    36353699#define __FUNCT__ "Penta::GetB_conduct"
    … …  
    36463710         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    36473711         */
    3648        
     3712
    36493713        int i;
    36503714        const int calculationdof=3;
    … …  
    36653729        }
    36663730}
    3667 
     3731/*}}}*/
     3732/*FUNCTION Penta GetB_artdiff {{{1*/
    36683733#undef __FUNCT__
    36693734#define __FUNCT__ "Penta::GetB_artdiff"
    … …  
    36973762        }
    36983763}
    3699 
     3764/*}}}*/
     3765/*FUNCTION Penta GetB_advec {{{1*/
    37003766#undef __FUNCT__
    37013767#define __FUNCT__ "Penta::GetB_advec"
    … …  
    37123778         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    37133779         */
    3714        
     3780
    37153781        int i;
    37163782        int calculationdof=3;
    … …  
    37313797        }
    37323798}
    3733 
     3799/*}}}*/
     3800/*FUNCTION Penta GetBprime_advec {{{1*/
    37343801#undef __FUNCT__
    37353802#define __FUNCT__ "Penta::GetBprime_advec"
    … …  
    37473814         * We assume B has been allocated already, of size: 3x(DOFPERGRID*numgrids)
    37483815         */
    3749        
     3816
    37503817        int i;
    37513818        const int calculationdof=3;
    … …  
    37663833        }
    37673834}
    3768 
     3835/*}}}*/
     3836/*FUNCTION Penta CreateKMatrixMelting {{{1*/
    37693837#undef __FUNCT__
    37703838#define __FUNCT__ "Penta::CreateKMatrixMelting"
    37713839void  Penta::CreateKMatrixMelting(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){
    3772        
     3840
    37733841        Tria* tria=NULL;
    3774        
     3842
    37753843        /*If on water, skip: */
    37763844        if(onwater)return;
    … …  
    37803848        }
    37813849        else{
    3782                
     3850
    37833851                tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
    37843852                tria->CreateKMatrixMelting(Kgg,inputs, analysis_type,sub_analysis_type);
    … …  
    37873855        }
    37883856}
    3789 
     3857/*}}}*/
     3858/*FUNCTION Penta CreatePVectorThermal {{{1*/
    37903859#undef __FUNCT__
    37913860#define __FUNCT__ "Penta::CreatePVectorThermal"
    … …  
    38053874        /*Grid data: */
    38063875        double        xyz_list[numgrids][3];
    3807        
     3876
    38083877        /* gaussian points: */
    38093878        int     num_area_gauss,igarea,igvert;
    … …  
    38183887        int     area_order=2;
    38193888        int       num_vert_gauss=3;
    3820        
     3889
    38213890        double dt;
    38223891        double vx_list[numgrids];
    … …  
    38343903        /* element parameters: */
    38353904        int    friction_type;
    3836        
     3905
    38373906        int    dofs[3]={0,1,2};
    38383907        int    dofs1[1]={0};
    … …  
    38723941        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    38733942        GetDofList(&doflist[0],&numberofdofspernode);
    3874        
     3943
    38753944        /*recovre material parameters: */
    38763945        rho_water=matpar->GetRhoWater();
    … …  
    38913960                if(!found)throw ErrorException(__FUNCT__," could not find temperature in inputs!");
    38923961        }
    3893        
     3962
    38943963        for(i=0;i<numgrids;i++){
    38953964                vx_list[i]=vxvyvz_list[i][0];
    … …  
    38993968
    39003969        /* Get gaussian points and weights. Penta is an extrusion of a Tria, we therefore
    3901            get tria gaussian points as well as segment gaussian points. For tria gaussian
    3902            points, order of integration is 2, because we need to integrate the product tB*D*B'
    3903            which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
    3904            points, same deal, which yields 3 gaussian points.: */
    3905        
     3970                get tria gaussian points as well as segment gaussian points. For tria gaussian
     3971                points, order of integration is 2, because we need to integrate the product tB*D*B'
     3972                which is a polynomial of degree 3 (see GaussTria for more details). For segment gaussian
     3973                points, same deal, which yields 3 gaussian points.: */
     3974
    39063975        /*Get gaussian points: */
    39073976        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);
    … …  
    39183987                        gauss_coord[2]=*(third_gauss_area_coord+igarea);
    39193988                        gauss_coord[3]=*(vert_gauss_coord+igvert);
    3920        
     3989
    39213990                        /*Compute strain rate and viscosity: */
    39223991                        GetStrainRateStokes(&epsilon[0],&vxvyvz_list[0][0],&xyz_list[0][0],gauss_coord);
    … …  
    39734042        extern int my_rank;
    39744043
    3975         cleanup_and_return:
     4044cleanup_and_return:
    39764045        xfree((void**)&first_gauss_area_coord);
    39774046        xfree((void**)&second_gauss_area_coord);
    … …  
    39824051
    39834052}
    3984 
     4053/*}}}*/
     4054/*FUNCTION Penta CreatePVectorMelting {{{1*/
    39854055#undef __FUNCT__
    39864056#define __FUNCT__ "Penta::CreatePVectorMelting"
    … …  
    39884058        return;
    39894059}
    3990 
     4060/*}}}*/
     4061/*FUNCTION Penta GetPhi {{{1*/
    39914062#undef __FUNCT__
    39924063#define __FUNCT__ "Penta::GetPhi"
    … …  
    40234094        *phi=2*pow(epsilon_eff,2.0)*viscosity;
    40244095}
    4025 
    4026 
     4096/*}}}*/
     4097/*FUNCTION Penta MassFlux {{{1*/
    40274098#undef __FUNCT__
    40284099#define __FUNCT__ "Penta::MassFlux"
    … …  
    40304101        throw ErrorException(__FUNCT__," not supported yet!");
    40314102}
     4103/*}}}*/

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

    r2705 r2711  
    19731973        *(J+NDOF2*1+1)=sqrt3/6.0*(2*y3-y1-y2);
    19741974}
    1975 
     1975/*}}}*/
     1976/*FUNCTION GetMatPar {{{1*/
    19761977void* Tria::GetMatPar(){
    19771978        return matpar;
    19781979}
    1979 
     1980/*}}}*/
     1981/*FUNCTION GetShelf {{{1*/
    19801982int   Tria::GetShelf(){
    19811983        return shelf;
    19821984}
    1983 
    1984                
     1985/*}}}*/
     1986/*FUNCTION GetNodes {{{1*/
    19851987void  Tria::GetNodes(void** vpnodes){
    19861988        int i;
    … …  
    19911993        }
    19921994}
    1993                
    1994 
     1995/*}}}*/
     1996/*FUNCTION GetOnBed {{{1*/
    19951997int Tria::GetOnBed(){
    19961998        return onbed;
    19971999}
    1998 
    1999 void          Tria::GetThicknessList(double* thickness_list){
     2000/*}}}*/
     2001/*FUNCTION GetThicknessList {{{1*/
     2002void Tria::GetThicknessList(double* thickness_list){
    20002003
    20012004        int i;
    20022005        for(i=0;i<3;i++)thickness_list[i]=h[i];
    20032006}
    2004 void          Tria::GetBedList(double* bed_list){
     2007/*}}}*/
     2008/*FUNCTION GetBedList {{{1*/
     2009void  Tria::GetBedList(double* bed_list){
    20052010       
    20062011        int i;
    … …  
    20082013
    20092014}
    2010  
    2011 
     2015/*}}}*/
     2016/*FUNCTION copy {{{1*/
    20122017Object* Tria::copy() {
    20132018       
Note: See TracChangeset for help on using the changeset viewer.