Changeset 2705


Ignore:
Timestamp:
12/07/09 16:24:28 (15 years ago)
Author:
Mathieu Morlighem
Message:

added folds in Tria.cpp

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2671 r2705  
    2424//#define _DEBUGGAUSS_
    2525
     26/*FUNCTION Tria constructor {{{1*/
    2627Tria::Tria(){
    2728        return;
    2829}
    29 
     30/*}}}*/
     31/*FUNCTION Tria destructor {{{1*/
    3032Tria::~Tria(){
    3133        return;
    3234}
    33 
     35/*}}}*/
     36/*FUNCTION Tria creation {{{1*/
    3437Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_numparid,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],
    3538                                double tria_accumulation[3],double tria_geothermalflux[3],int tria_friction_type,double tria_p,double tria_q,int tria_shelf, bool tria_onwater){
    … …  
    6871        return;
    6972}
    70 
     73/*}}}*/
     74/*FUNCTION Echo {{{1*/
    7175#undef __FUNCT__
    7276#define __FUNCT__ "Tria::Echo"
    … …  
    104108        return;
    105109}
     110/*}}}*/
     111/*FUNCTION DeepEcho{{{1*/
    106112#undef __FUNCT__
    107113#define __FUNCT__ "Tria::DeepEcho"
    … …  
    139145        return;
    140146}
     147/*}}}*/
     148/*FUNCTION Marshall {{{1*/
    141149void  Tria::Marshall(char** pmarshalled_dataset){
    142150
    … …  
    184192        return;
    185193}
    186                
     194/*}}}*/
     195/*FUNCTION MarshallSize {{{1*/
    187196int   Tria::MarshallSize(){
    188197        return sizeof(id)
    … …  
    214223                +sizeof(int); //sizeof(int) for enum type
    215224}
    216 
     225/*}}}*/
     226/*FUNCTION GetName {{{1*/
    217227char* Tria::GetName(void){
    218228        return "tria";
    219229}
    220 
     230/*}}}*/
     231/*FUNCTION Demarshall {{{1*/
    221232void  Tria::Demarshall(char** pmarshalled_dataset){
    222233
    … …  
    267278        return;
    268279}
     280/*}}}*/
     281/*FUNCTION Enum {{{1*/
    269282int Tria::Enum(void){
    270283
    … …  
    272285
    273286}
     287/*}}}*/
     288/*FUNCTION GetId {{{1*/
    274289int    Tria::GetId(){ return id; }
    275 
     290/*}}}*/
     291/*FUNCTION MyRank {{{1*/
    276292int    Tria::MyRank(void){
    277293        extern int my_rank;
    278294        return my_rank;
    279295}
    280 
    281 
     296/*}}}*/
     297/*FUNCTION Configure {{{1*/
    282298#undef __FUNCT__
    283299#define __FUNCT__ "Tria::Configure"
    … …  
    308324
    309325}
    310 
     326/*}}}*/
     327/*FUNCTION CreateKMatrix {{{1*/
    311328#undef __FUNCT__
    312329#define __FUNCT__ "Tria::CreateKMatrix"
    … …  
    343360
    344361}
    345 
    346 
     362/*}}}*/
     363/*FUNCTION CreateKMatrixDiagnosticHoriz {{{1*/
    347364#undef __FUNCT__
    348365#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz"
    … …  
    562579
    563580}
     581/*}}}*/
     582/*FUNCTION CreateKMatrixPrognostic {{{1*/
    564583#undef __FUNCT__
    565584#define __FUNCT__ "Tria::CreateKMatrixPrognostic"
    … …  
    780799
    781800}
    782 
    783 #undef __FUNCT__
    784 #define __FUNCT__ "Tria::CreatePVectorPrognostic"
    785 void  Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
     801/*}}}*/
     802/*FUNCTION CreateKMatrixDiagnosticHorizFriction {{{1*/
     803#undef __FUNCT__
     804#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction"
     805void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    786806
    787807
    … …  
    791811        /* node data: */
    792812        const int    numgrids=3;
    793         const int    NDOF1=1;
    794         const int    numdof=NDOF1*numgrids;
     813        const int    numdof=2*numgrids;
    795814        double       xyz_list[numgrids][3];
    796815        int          doflist[numdof];
    … …  
    806825        double  gauss_l1l2l3[3];
    807826
     827        /* matrices: */
     828        double L[2][numdof];
     829        double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
     830        double DL_scalar;
     831
     832        /* local element matrices: */
     833        double Ke_gg[numdof][numdof]; //local element stiffness matrix
     834        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
     835       
     836        double Jdet;
     837       
     838        /*slope: */
     839        double  slope[2]={0.0,0.0};
     840        double  slope_magnitude;
     841
     842        /*input parameters for structural analysis (diagnostic): */
     843        double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
     844        int     dofs[2]={0,1};
     845
     846        /*friction: */
     847        double alpha2_list[numgrids]={0.0,0.0,0.0};
     848        double alpha2;
     849
     850        double MAXSLOPE=.06; // 6 %
     851        double MOUNTAINKEXPONENT=10;
     852
     853        ParameterInputs* inputs=NULL;
     854
     855        /*recover pointers: */
     856        inputs=(ParameterInputs*)vinputs;
     857       
     858        /* Get node coordinates and dof list: */
     859        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     860        GetDofList(&doflist[0],&numberofdofspernode);
     861
     862        /* Set Ke_gg to 0: */
     863        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
     864
     865        if (shelf){
     866                /*no friction, do nothing*/
     867                return;
     868        }
     869
     870        if (friction_type!=2)throw ErrorException(__FUNCT__," non-viscous friction not supported yet!");
     871
     872        /*recover extra inputs from users, at current convergence iteration: */
     873        inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
     874
     875        /*Build alpha2_list used by drag stiffness matrix*/
     876        Friction* friction=NewFriction();
     877       
     878        /*Initialize all fields: */
     879        friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
     880        strcpy(friction->element_type,"2d");
     881       
     882        friction->gravity=matpar->GetG();
     883        friction->rho_ice=matpar->GetRhoIce();
     884        friction->rho_water=matpar->GetRhoWater();
     885        friction->K=&k[0];
     886        friction->bed=&b[0];
     887        friction->thickness=&h[0];
     888        friction->velocities=&vxvy_list[0][0];
     889        friction->p=p;
     890        friction->q=q;
     891
     892        /*Compute alpha2_list: */
     893        FrictionGetAlpha2(&alpha2_list[0],friction);
     894
     895        /*Erase friction object: */
     896        DeleteFriction(&friction);
     897
     898        #ifdef _DEBUGELEMENTS_
     899        if(my_rank==RANK && id==ELID){
     900                printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
     901        }
     902        #endif
     903
     904        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     905        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     906
     907        #ifdef _DEBUGELEMENTS_
     908        if(my_rank==RANK && id==ELID){
     909                printf("   gaussian points: \n");
     910                for(i=0;i<num_gauss;i++){
     911                        printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
     912                }
     913        }
     914        #endif
     915
     916        /* Start  looping on the number of gaussian points: */
     917        for (ig=0; ig<num_gauss; ig++){
     918                /*Pick up the gaussian point: */
     919                gauss_weight=*(gauss_weights+ig);
     920                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     921                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     922                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     923
     924
     925                // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
     926                //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
     927                GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
     928                slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
     929
     930                if (slope_magnitude>MAXSLOPE){
     931                        alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT);
     932                        alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT);
     933                        alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT);
     934                }
     935
     936                /* Get Jacobian determinant: */
     937                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     938
     939                /*Get L matrix: */
     940                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
     941
     942                /*Now, take care of the basal friction if there is any: */
     943                GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
     944
     945                DL_scalar=alpha2*gauss_weight*Jdet;
     946                for (i=0;i<2;i++){
     947                        DL[i][i]=DL_scalar;
     948                }
     949               
     950                /*  Do the triple producte tL*D*L: */
     951                TripleMultiply( &L[0][0],2,numdof,1,
     952                                        &DL[0][0],2,2,0,
     953                                        &L[0][0],2,numdof,0,
     954                                        &Ke_gg_gaussian[0][0],0);
     955
     956                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     957
     958        } // for (ig=0; ig<num_gauss; ig++)
     959
     960        /*Add Ke_gg to global matrix Kgg: */
     961        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     962
     963        cleanup_and_return:
     964        xfree((void**)&first_gauss_area_coord);
     965        xfree((void**)&second_gauss_area_coord);
     966        xfree((void**)&third_gauss_area_coord);
     967        xfree((void**)&gauss_weights);
     968
     969}       
     970/*}}}*/
     971/*FUNCTION CreateKMatrixSlopeCompute {{{1*/
     972#undef __FUNCT__
     973#define __FUNCT__ "Tria::CreateKMatrixSlopeCompute"
     974
     975void  Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
     976
     977        /* local declarations */
     978        int             i,j;
     979
     980        /* node data: */
     981        const int    numgrids=3;
     982        const int    NDOF1=1;
     983        const int    numdof=NDOF1*numgrids;
     984        double       xyz_list[numgrids][3];
     985        int          doflist[numdof];
     986        int          numberofdofspernode;
     987       
     988        /* gaussian points: */
     989        int     num_gauss,ig;
     990        double* first_gauss_area_coord  =  NULL;
     991        double* second_gauss_area_coord =  NULL;
     992        double* third_gauss_area_coord  =  NULL;
     993        double* gauss_weights           =  NULL;
     994        double  gauss_weight;
     995        double  gauss_l1l2l3[3];
     996
     997        /* matrices: */
     998        double L[1][3];
     999        double DL_scalar;
     1000
     1001        /* local element matrices: */
     1002        double Ke_gg[numdof][numdof]; //local element stiffness matrix
     1003        double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
     1004       
     1005        double Jdet;
     1006       
     1007        /* Get node coordinates and dof list: */
     1008        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     1009        GetDofList(&doflist[0],&numberofdofspernode);
     1010
     1011        /* Set Ke_gg to 0: */
     1012        for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
     1013
     1014        /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
     1015        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
     1016
     1017        /* Start  looping on the number of gaussian points: */
     1018        for (ig=0; ig<num_gauss; ig++){
     1019                /*Pick up the gaussian point: */
     1020                gauss_weight=*(gauss_weights+ig);
     1021                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1022                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1023                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1024
     1025               
     1026                /*Get L matrix: */
     1027                GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
     1028
     1029                /* Get Jacobian determinant: */
     1030                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1031               
     1032                DL_scalar=gauss_weight*Jdet;
     1033
     1034                /*  Do the triple producte tL*D*L: */
     1035                TripleMultiply( &L[0][0],1,3,1,
     1036                                        &DL_scalar,1,1,0,
     1037                                        &L[0][0],1,3,0,
     1038                                        &Ke_gg_gaussian[0][0],0);
     1039
     1040                /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
     1041                for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
     1042        } //for (ig=0; ig<num_gauss; ig++
     1043
     1044        /*Add Ke_gg to global matrix Kgg: */
     1045        MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     1046               
     1047        cleanup_and_return:
     1048        xfree((void**)&first_gauss_area_coord);
     1049        xfree((void**)&second_gauss_area_coord);
     1050        xfree((void**)&third_gauss_area_coord);
     1051        xfree((void**)&gauss_weights);
     1052}
     1053/*}}}*/
     1054/*FUNCTION CreatePVector {{{1*/
     1055#undef __FUNCT__
     1056#define __FUNCT__ "Tria::CreatePVector"
     1057void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
     1058       
     1059        /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
     1060        if (analysis_type==ControlAnalysisEnum()){
     1061               
     1062                CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     1063       
     1064        }
     1065        else if (analysis_type==DiagnosticAnalysisEnum()){
     1066                if (sub_analysis_type==HorizAnalysisEnum()){
     1067               
     1068                        CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
     1069               
     1070                }
     1071                else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
     1072        }
     1073        else if (analysis_type==SlopeComputeAnalysisEnum()){
     1074               
     1075                CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
     1076        }
     1077        else if (analysis_type==PrognosticAnalysisEnum()){
     1078
     1079                CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
     1080        }
     1081        else{
     1082                throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
     1083        }
     1084
     1085}
     1086/*}}}*/
     1087/*FUNCTION CreatePVectorDiagnosticHoriz {{{1*/
     1088#undef __FUNCT__
     1089#define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"
     1090void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
     1091
     1092        int             i,j;
     1093
     1094        /* node data: */
     1095        const int    numgrids=3;
     1096        const int    numdof=2*numgrids;
     1097        const int    NDOF2=2;
     1098        double       xyz_list[numgrids][3];
     1099        int          doflist[numdof];
     1100        int          numberofdofspernode;
     1101       
     1102        /* parameters: */
     1103        double  plastic_stress;
     1104        double  slope[NDOF2];
     1105        double  driving_stress_baseline;
     1106
     1107        /* gaussian points: */
     1108        int     num_gauss,ig;
     1109        double* first_gauss_area_coord  =  NULL;
     1110        double* second_gauss_area_coord =  NULL;
     1111        double* third_gauss_area_coord  =  NULL;
     1112        double* gauss_weights           =  NULL;
     1113        double  gauss_weight;
     1114        double  gauss_l1l2l3[3];
     1115
     1116        /* Jacobian: */
     1117        double Jdet;
     1118
     1119        /*nodal functions: */
     1120        double l1l2l3[3];
     1121
     1122        /*element vector at the gaussian points: */
     1123        double  pe_g[numdof];
     1124        double  pe_g_gaussian[numdof];
     1125
     1126        /*input parameters for structural analysis (diagnostic): */
     1127        double  thickness;
     1128
     1129        ParameterInputs* inputs=NULL;
     1130
     1131        /*First, if we are on water, return empty vector: */
     1132        if(onwater)return;
     1133
     1134        /*recover pointers: */
     1135        inputs=(ParameterInputs*)vinputs;
     1136
     1137        /* Get node coordinates and dof list: */
     1138        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
     1139        GetDofList(&doflist[0],&numberofdofspernode);
     1140
     1141        /* Set pe_g to 0: */
     1142        for(i=0;i<numdof;i++) pe_g[i]=0.0;
     1143
     1144
     1145        #ifdef _DEBUGELEMENTS_
     1146        if(my_rank==RANK && id==ELID){
     1147                printf("gravity %g\n",matpar->GetG());
     1148                printf("rho_ice %g\n",matpar->GetRhoIce());
     1149                printf("thickness [%g,%g,%g]\n",h[0],h[1],h[2]);
     1150                printf("surface[%g,%g,%g]\n",s[0],s[1],s[2]);
     1151                printf("bed[%g,%g,%g]\n",b[0],b[1],b[2]);
     1152                printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]);
     1153        }
     1154        #endif
     1155
     1156
     1157        /* Get gaussian points and weights: */
     1158        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
     1159
     1160        #ifdef _DEBUGELEMENTS_
     1161        if(my_rank==RANK && id==ELID){
     1162                printf("   gaussian points: \n");
     1163                for(i=0;i<num_gauss;i++){
     1164                        printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
     1165                }
     1166        }
     1167        #endif
     1168
     1169
     1170
     1171        /* Start  looping on the number of gaussian points: */
     1172        for (ig=0; ig<num_gauss; ig++){
     1173                /*Pick up the gaussian point: */
     1174                gauss_weight=*(gauss_weights+ig);
     1175                gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
     1176                gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
     1177                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
     1178
     1179                /*Compute thickness at gaussian point: */
     1180                GetParameterValue(&thickness, &h[0],gauss_l1l2l3);
     1181       
     1182                GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
     1183               
     1184                /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the
     1185                 * element itself: */
     1186                if(friction_type==1){
     1187                        GetParameterValue(&plastic_stress, &k[0],gauss_l1l2l3);
     1188                }
     1189
     1190                /* Get Jacobian determinant: */
     1191                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
     1192               
     1193                 /*Get nodal functions: */
     1194                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     1195
     1196                /*Compute driving stress: */
     1197                driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
     1198
     1199
     1200                #ifdef _DEBUGELEMENTS_
     1201                if(my_rank==RANK && id==ELID){
     1202                        printf("      gaussian %i\n",ig);
     1203                        printf("      thickness %g\n",thickness);
     1204                        printf("      slope(%g,%g)\n",slope[0],slope[1]);
     1205                        printf("      Jdet %g\n",Jdet);
     1206                        printf("      gaussweigth %g\n",gauss_weight);
     1207                        printf("      l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]);
     1208                        if(friction_type==1)printf("      plastic_stress(%g)\n",plastic_stress);
     1209                }
     1210                #endif
     1211
     1212                /*Build pe_g_gaussian vector: */
     1213                if(friction_type==1){
     1214                        for (i=0;i<numgrids;i++){
     1215                                for (j=0;j<NDOF2;j++){
     1216                                        pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss_weight*l1l2l3[i];
     1217                                }
     1218                        }
     1219                }
     1220                else {
     1221                        for (i=0;i<numgrids;i++){
     1222                                for (j=0;j<NDOF2;j++){
     1223                                        pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3[i];
     1224                                }
     1225                        }
     1226                }
     1227
     1228                /*Add pe_g_gaussian vector to pe_g: */
     1229                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
     1230
     1231        } //for (ig=0; ig<num_gauss; ig++)
     1232
     1233        #ifdef _DEBUGELEMENTS_
     1234        if(my_rank==RANK && id==ELID){
     1235                printf("      pe_g->terms\n",ig);
     1236                for( i=0; i<pe_g->nrows; i++){
     1237                        printf("%g ",*(pe_g->terms+i));
     1238                }
     1239                printf("\n");
     1240                printf("      pe_g->row_indices\n",ig);
     1241                for( i=0; i<pe_g->nrows; i++){
     1242                        printf("%i ",*(pe_g->row_indices+i));
     1243                }
     1244        }
     1245        #endif
     1246
     1247        /*Add pe_g to global vector pg: */
     1248        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
     1249
     1250        cleanup_and_return:
     1251        xfree((void**)&first_gauss_area_coord);
     1252        xfree((void**)&second_gauss_area_coord);
     1253        xfree((void**)&third_gauss_area_coord);
     1254        xfree((void**)&gauss_weights);
     1255
     1256}
     1257/*}}}*/
     1258/*FUNCTION CreatePVectorPrognostic {{{1*/
     1259#undef __FUNCT__
     1260#define __FUNCT__ "Tria::CreatePVectorPrognostic"
     1261void  Tria::CreatePVectorPrognostic(Vec pg ,void* vinputs,int analysis_type,int sub_analysis_type){
     1262
     1263
     1264        /* local declarations */
     1265        int             i,j;
     1266
     1267        /* node data: */
     1268        const int    numgrids=3;
     1269        const int    NDOF1=1;
     1270        const int    numdof=NDOF1*numgrids;
     1271        double       xyz_list[numgrids][3];
     1272        int          doflist[numdof];
     1273        int          numberofdofspernode;
     1274
     1275        /* gaussian points: */
     1276        int     num_gauss,ig;
     1277        double* first_gauss_area_coord  =  NULL;
     1278        double* second_gauss_area_coord =  NULL;
     1279        double* third_gauss_area_coord  =  NULL;
     1280        double* gauss_weights           =  NULL;
     1281        double  gauss_weight;
     1282        double  gauss_l1l2l3[3];
     1283
    8081284        /* matrix */
    8091285        double pe_g[numgrids]={0.0};
    8101286        double L[numgrids];
    8111287        double Jdettria;
    812        
     1288
    8131289        /*input parameters for structural analysis (diagnostic): */
    8141290        double  accumulation_list[numgrids]={0.0};
    … …  
    8621338                GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3);
    8631339                GetParameterValue(&thickness_g, &thickness_list[0],gauss_l1l2l3);
    864                
     1340
    8651341                /* Add value into pe_g: */
    8661342                for( i=0; i<numdof; i++) pe_g[i]+=Jdettria*gauss_weight*(thickness_g+dt*(accumulation_g-melting_g))*L[i];
    867                
     1343
    8681344        } // for (ig=0; ig<num_gauss; ig++)
    8691345
    … …  
    8711347        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    8721348
    873         cleanup_and_return:
     1349cleanup_and_return:
    8741350        xfree((void**)&first_gauss_area_coord);
    8751351        xfree((void**)&second_gauss_area_coord);
    … …  
    8781354
    8791355}
    880 
    881 
    882 #undef __FUNCT__
    883 #define __FUNCT__ "Tria::CreateKMatrixDiagnosticHorizFriction"
    884 void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    885 
    886 
    887         /* local declarations */
     1356/*}}}*/
     1357/*FUNCTION CreatePVectorSlopeCompute {{{1*/
     1358#undef __FUNCT__
     1359#define __FUNCT__ "Tria::CreatePVectorSlopeCompute"
     1360
     1361void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
     1362
    8881363        int             i,j;
    8891364
    8901365        /* node data: */
    8911366        const int    numgrids=3;
    892         const int    numdof=2*numgrids;
     1367        const int    NDOF1=1;
     1368        const int    numdof=NDOF1*numgrids;
    8931369        double       xyz_list[numgrids][3];
    8941370        int          doflist[numdof];
    … …  
    9041380        double  gauss_l1l2l3[3];
    9051381
    906         /* matrices: */
    907         double L[2][numdof];
    908         double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
    909         double DL_scalar;
    910 
    911         /* local element matrices: */
    912         double Ke_gg[numdof][numdof]; //local element stiffness matrix
    913         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
    914        
     1382        /* Jacobian: */
    9151383        double Jdet;
    916        
    917         /*slope: */
    918         double  slope[2]={0.0,0.0};
    919         double  slope_magnitude;
    920 
    921         /*input parameters for structural analysis (diagnostic): */
    922         double  vxvy_list[numgrids][2]={{0,0},{0,0},{0,0}};
    923         int     dofs[2]={0,1};
    924 
    925         /*friction: */
    926         double alpha2_list[numgrids]={0.0,0.0,0.0};
    927         double alpha2;
    928 
    929         double MAXSLOPE=.06; // 6 %
    930         double MOUNTAINKEXPONENT=10;
    931 
    932         ParameterInputs* inputs=NULL;
    933 
    934         /*recover pointers: */
    935         inputs=(ParameterInputs*)vinputs;
    936        
     1384
     1385        /*nodal functions: */
     1386        double l1l2l3[3];
     1387
     1388        /*element vector at the gaussian points: */
     1389        double  pe_g[numdof];
     1390        double  pe_g_gaussian[numdof];
     1391        double  param[3];
     1392        double  slope[2];
     1393
    9371394        /* Get node coordinates and dof list: */
    9381395        GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    9391396        GetDofList(&doflist[0],&numberofdofspernode);
    9401397
    941         /* Set Ke_gg to 0: */
    942         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    943 
    944         if (shelf){
    945                 /*no friction, do nothing*/
    946                 return;
    947         }
    948 
    949         if (friction_type!=2)throw ErrorException(__FUNCT__," non-viscous friction not supported yet!");
    950 
    951         /*recover extra inputs from users, at current convergence iteration: */
    952         inputs->Recover("velocity",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes);
    953 
    954         /*Build alpha2_list used by drag stiffness matrix*/
    955         Friction* friction=NewFriction();
    956        
    957         /*Initialize all fields: */
    958         friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
    959         strcpy(friction->element_type,"2d");
    960        
    961         friction->gravity=matpar->GetG();
    962         friction->rho_ice=matpar->GetRhoIce();
    963         friction->rho_water=matpar->GetRhoWater();
    964         friction->K=&k[0];
    965         friction->bed=&b[0];
    966         friction->thickness=&h[0];
    967         friction->velocities=&vxvy_list[0][0];
    968         friction->p=p;
    969         friction->q=q;
    970 
    971         /*Compute alpha2_list: */
    972         FrictionGetAlpha2(&alpha2_list[0],friction);
    973 
    974         /*Erase friction object: */
    975         DeleteFriction(&friction);
    976 
    977         #ifdef _DEBUGELEMENTS_
    978         if(my_rank==RANK && id==ELID){
    979                 printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
    980         }
    981         #endif
    982 
    983         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    984         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    985 
    986         #ifdef _DEBUGELEMENTS_
    987         if(my_rank==RANK && id==ELID){
    988                 printf("   gaussian points: \n");
    989                 for(i=0;i<num_gauss;i++){
    990                         printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
    991                 }
    992         }
    993         #endif
     1398        /* Set pe_g to 0: */
     1399        for(i=0;i<numdof;i++) pe_g[i]=0.0;
     1400
     1401        if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==SurfaceYAnalysisEnum())){
     1402                for(i=0;i<numdof;i++) param[i]=s[i];
     1403        }
     1404        if ( (sub_analysis_type==BedXAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
     1405                for(i=0;i<numdof;i++) param[i]=b[i];
     1406        }
     1407
     1408        /* Get gaussian points and weights: */
     1409        GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
     1410
    9941411
    9951412        /* Start  looping on the number of gaussian points: */
    … …  
    10011418                gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    10021419
    1003 
    1004                 // If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case,
    1005                 //velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
    1006                 GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
    1007                 slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
    1008 
    1009                 if (slope_magnitude>MAXSLOPE){
    1010                         alpha2_list[0]=pow((double)10,MOUNTAINKEXPONENT);
    1011                         alpha2_list[1]=pow((double)10,MOUNTAINKEXPONENT);
    1012                         alpha2_list[2]=pow((double)10,MOUNTAINKEXPONENT);
    1013                 }
    1014 
     1420                GetParameterDerivativeValue(&slope[0], &param[0],&xyz_list[0][0], gauss_l1l2l3);
     1421               
    10151422                /* Get Jacobian determinant: */
    10161423                GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1017 
    1018                 /*Get L matrix: */
    1019                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
    1020 
    1021                 /*Now, take care of the basal friction if there is any: */
    1022                 GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);
    1023 
    1024                 DL_scalar=alpha2*gauss_weight*Jdet;
    1025                 for (i=0;i<2;i++){
    1026                         DL[i][i]=DL_scalar;
    1027                 }
    1028                
    1029                 /*  Do the triple producte tL*D*L: */
    1030                 TripleMultiply( &L[0][0],2,numdof,1,
    1031                                         &DL[0][0],2,2,0,
    1032                                         &L[0][0],2,numdof,0,
    1033                                         &Ke_gg_gaussian[0][0],0);
    1034 
    1035                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    1036 
    1037         } // for (ig=0; ig<num_gauss; ig++)
    1038 
    1039         /*Add Ke_gg to global matrix Kgg: */
    1040         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
     1424               
     1425                 /*Get nodal functions: */
     1426                GetNodalFunctions(l1l2l3, gauss_l1l2l3);
     1427
     1428                /*Build pe_g_gaussian vector: */
     1429                if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==BedXAnalysisEnum())){
     1430                        for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
     1431                }
     1432                if ( (sub_analysis_type==SurfaceYAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
     1433                        for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
     1434                }
     1435
     1436                /*Add pe_g_gaussian vector to pe_g: */
     1437                for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
     1438
     1439        } //for (ig=0; ig<num_gauss; ig++)
     1440
     1441        /*Add pe_g to global vector pg: */
     1442        VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    10411443
    10421444        cleanup_and_return:
    … …  
    10461448        xfree((void**)&gauss_weights);
    10471449
    1048 }       
    1049 
    1050 #undef __FUNCT__
    1051 #define __FUNCT__ "Tria::CreateKMatrixSlopeCompute"
    1052 
    1053 void  Tria::CreateKMatrixSlopeCompute(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type){
    1054 
    1055         /* local declarations */
    1056         int             i,j;
    1057 
    1058         /* node data: */
    1059         const int    numgrids=3;
    1060         const int    NDOF1=1;
    1061         const int    numdof=NDOF1*numgrids;
    1062         double       xyz_list[numgrids][3];
    1063         int          doflist[numdof];
    1064         int          numberofdofspernode;
    1065        
    1066         /* gaussian points: */
    1067         int     num_gauss,ig;
    1068         double* first_gauss_area_coord  =  NULL;
    1069         double* second_gauss_area_coord =  NULL;
    1070         double* third_gauss_area_coord  =  NULL;
    1071         double* gauss_weights           =  NULL;
    1072         double  gauss_weight;
    1073         double  gauss_l1l2l3[3];
    1074 
    1075         /* matrices: */
    1076         double L[1][3];
    1077         double DL_scalar;
    1078 
    1079         /* local element matrices: */
    1080         double Ke_gg[numdof][numdof]; //local element stiffness matrix
    1081         double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
    1082        
    1083         double Jdet;
    1084        
    1085         /* Get node coordinates and dof list: */
    1086         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    1087         GetDofList(&doflist[0],&numberofdofspernode);
    1088 
    1089         /* Set Ke_gg to 0: */
    1090         for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke_gg[i][j]=0.0;
    1091 
    1092         /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
    1093         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
    1094 
    1095         /* Start  looping on the number of gaussian points: */
    1096         for (ig=0; ig<num_gauss; ig++){
    1097                 /*Pick up the gaussian point: */
    1098                 gauss_weight=*(gauss_weights+ig);
    1099                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1100                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1101                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1102 
    1103                
    1104                 /*Get L matrix: */
    1105                 GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,NDOF1);
    1106 
    1107                 /* Get Jacobian determinant: */
    1108                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1109                
    1110                 DL_scalar=gauss_weight*Jdet;
    1111 
    1112                 /*  Do the triple producte tL*D*L: */
    1113                 TripleMultiply( &L[0][0],1,3,1,
    1114                                         &DL_scalar,1,1,0,
    1115                                         &L[0][0],1,3,0,
    1116                                         &Ke_gg_gaussian[0][0],0);
    1117 
    1118                 /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
    1119                 for( i=0; i<numdof; i++) for(j=0;j<numdof;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
    1120         } //for (ig=0; ig<num_gauss; ig++
    1121 
    1122         /*Add Ke_gg to global matrix Kgg: */
    1123         MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
    1124                
    1125         cleanup_and_return:
    1126         xfree((void**)&first_gauss_area_coord);
    1127         xfree((void**)&second_gauss_area_coord);
    1128         xfree((void**)&third_gauss_area_coord);
    1129         xfree((void**)&gauss_weights);
    1130 }
    1131 
    1132 #undef __FUNCT__
    1133 #define __FUNCT__ "Tria::CreatePVector"
    1134 void  Tria::CreatePVector(Vec pg,void* inputs,int analysis_type,int sub_analysis_type){
    1135        
    1136         /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    1137         if (analysis_type==ControlAnalysisEnum()){
    1138                
    1139                 CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
    1140        
    1141         }
    1142         else if (analysis_type==DiagnosticAnalysisEnum()){
    1143                 if (sub_analysis_type==HorizAnalysisEnum()){
    1144                
    1145                         CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type,sub_analysis_type);
    1146                
    1147                 }
    1148                 else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet"));
    1149         }
    1150         else if (analysis_type==SlopeComputeAnalysisEnum()){
    1151                
    1152                 CreatePVectorSlopeCompute( pg,inputs,analysis_type,sub_analysis_type);
    1153         }
    1154         else if (analysis_type==PrognosticAnalysisEnum()){
    1155 
    1156                 CreatePVectorPrognostic( pg,inputs,analysis_type,sub_analysis_type);
    1157         }
    1158         else{
    1159                 throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
    1160         }
    1161 
    1162 }
    1163 
    1164 #undef __FUNCT__
    1165 #define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"
    1166 void Tria::CreatePVectorDiagnosticHoriz( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
    1167 
    1168         int             i,j;
    1169 
    1170         /* node data: */
    1171         const int    numgrids=3;
    1172         const int    numdof=2*numgrids;
    1173         const int    NDOF2=2;
    1174         double       xyz_list[numgrids][3];
    1175         int          doflist[numdof];
    1176         int          numberofdofspernode;
    1177        
    1178         /* parameters: */
    1179         double  plastic_stress;
    1180         double  slope[NDOF2];
    1181         double  driving_stress_baseline;
    1182 
    1183         /* gaussian points: */
    1184         int     num_gauss,ig;
    1185         double* first_gauss_area_coord  =  NULL;
    1186         double* second_gauss_area_coord =  NULL;
    1187         double* third_gauss_area_coord  =  NULL;
    1188         double* gauss_weights           =  NULL;
    1189         double  gauss_weight;
    1190         double  gauss_l1l2l3[3];
    1191 
    1192         /* Jacobian: */
    1193         double Jdet;
    1194 
    1195         /*nodal functions: */
    1196         double l1l2l3[3];
    1197 
    1198         /*element vector at the gaussian points: */
    1199         double  pe_g[numdof];
    1200         double  pe_g_gaussian[numdof];
    1201 
    1202         /*input parameters for structural analysis (diagnostic): */
    1203         double  thickness;
    1204 
    1205         ParameterInputs* inputs=NULL;
    1206 
    1207         /*First, if we are on water, return empty vector: */
    1208         if(onwater)return;
    1209 
    1210         /*recover pointers: */
    1211         inputs=(ParameterInputs*)vinputs;
    1212 
    1213         /* Get node coordinates and dof list: */
    1214         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    1215         GetDofList(&doflist[0],&numberofdofspernode);
    1216 
    1217         /* Set pe_g to 0: */
    1218         for(i=0;i<numdof;i++) pe_g[i]=0.0;
    1219 
    1220 
    1221         #ifdef _DEBUGELEMENTS_
    1222         if(my_rank==RANK && id==ELID){
    1223                 printf("gravity %g\n",matpar->GetG());
    1224                 printf("rho_ice %g\n",matpar->GetRhoIce());
    1225                 printf("thickness [%g,%g,%g]\n",h[0],h[1],h[2]);
    1226                 printf("surface[%g,%g,%g]\n",s[0],s[1],s[2]);
    1227                 printf("bed[%g,%g,%g]\n",b[0],b[1],b[2]);
    1228                 printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]);
    1229         }
    1230         #endif
    1231 
    1232 
    1233         /* Get gaussian points and weights: */
    1234         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
    1235 
    1236         #ifdef _DEBUGELEMENTS_
    1237         if(my_rank==RANK && id==ELID){
    1238                 printf("   gaussian points: \n");
    1239                 for(i=0;i<num_gauss;i++){
    1240                         printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
    1241                 }
    1242         }
    1243         #endif
    1244 
    1245 
    1246 
    1247         /* Start  looping on the number of gaussian points: */
    1248         for (ig=0; ig<num_gauss; ig++){
    1249                 /*Pick up the gaussian point: */
    1250                 gauss_weight=*(gauss_weights+ig);
    1251                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1252                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1253                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1254 
    1255                 /*Compute thickness at gaussian point: */
    1256                 GetParameterValue(&thickness, &h[0],gauss_l1l2l3);
    1257        
    1258                 GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
    1259                
    1260                 /*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the
    1261                  * element itself: */
    1262                 if(friction_type==1){
    1263                         GetParameterValue(&plastic_stress, &k[0],gauss_l1l2l3);
    1264                 }
    1265 
    1266                 /* Get Jacobian determinant: */
    1267                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1268                
    1269                  /*Get nodal functions: */
    1270                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    1271 
    1272                 /*Compute driving stress: */
    1273                 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
    1274 
    1275 
    1276                 #ifdef _DEBUGELEMENTS_
    1277                 if(my_rank==RANK && id==ELID){
    1278                         printf("      gaussian %i\n",ig);
    1279                         printf("      thickness %g\n",thickness);
    1280                         printf("      slope(%g,%g)\n",slope[0],slope[1]);
    1281                         printf("      Jdet %g\n",Jdet);
    1282                         printf("      gaussweigth %g\n",gauss_weight);
    1283                         printf("      l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]);
    1284                         if(friction_type==1)printf("      plastic_stress(%g)\n",plastic_stress);
    1285                 }
    1286                 #endif
    1287 
    1288                 /*Build pe_g_gaussian vector: */
    1289                 if(friction_type==1){
    1290                         for (i=0;i<numgrids;i++){
    1291                                 for (j=0;j<NDOF2;j++){
    1292                                         pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss_weight*l1l2l3[i];
    1293                                 }
    1294                         }
    1295                 }
    1296                 else {
    1297                         for (i=0;i<numgrids;i++){
    1298                                 for (j=0;j<NDOF2;j++){
    1299                                         pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3[i];
    1300                                 }
    1301                         }
    1302                 }
    1303 
    1304                 /*Add pe_g_gaussian vector to pe_g: */
    1305                 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    1306 
    1307         } //for (ig=0; ig<num_gauss; ig++)
    1308 
    1309         #ifdef _DEBUGELEMENTS_
    1310         if(my_rank==RANK && id==ELID){
    1311                 printf("      pe_g->terms\n",ig);
    1312                 for( i=0; i<pe_g->nrows; i++){
    1313                         printf("%g ",*(pe_g->terms+i));
    1314                 }
    1315                 printf("\n");
    1316                 printf("      pe_g->row_indices\n",ig);
    1317                 for( i=0; i<pe_g->nrows; i++){
    1318                         printf("%i ",*(pe_g->row_indices+i));
    1319                 }
    1320         }
    1321         #endif
    1322 
    1323         /*Add pe_g to global vector pg: */
    1324         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    1325 
    1326         cleanup_and_return:
    1327         xfree((void**)&first_gauss_area_coord);
    1328         xfree((void**)&second_gauss_area_coord);
    1329         xfree((void**)&third_gauss_area_coord);
    1330         xfree((void**)&gauss_weights);
    1331 
    1332 }
    1333 
    1334 #undef __FUNCT__
    1335 #define __FUNCT__ "Tria::CreatePVectorSlopeCompute"
    1336 
    1337 void Tria::CreatePVectorSlopeCompute( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type){
    1338 
    1339         int             i,j;
    1340 
    1341         /* node data: */
    1342         const int    numgrids=3;
    1343         const int    NDOF1=1;
    1344         const int    numdof=NDOF1*numgrids;
    1345         double       xyz_list[numgrids][3];
    1346         int          doflist[numdof];
    1347         int          numberofdofspernode;
    1348        
    1349         /* gaussian points: */
    1350         int     num_gauss,ig;
    1351         double* first_gauss_area_coord  =  NULL;
    1352         double* second_gauss_area_coord =  NULL;
    1353         double* third_gauss_area_coord  =  NULL;
    1354         double* gauss_weights           =  NULL;
    1355         double  gauss_weight;
    1356         double  gauss_l1l2l3[3];
    1357 
    1358         /* Jacobian: */
    1359         double Jdet;
    1360 
    1361         /*nodal functions: */
    1362         double l1l2l3[3];
    1363 
    1364         /*element vector at the gaussian points: */
    1365         double  pe_g[numdof];
    1366         double  pe_g_gaussian[numdof];
    1367         double  param[3];
    1368         double  slope[2];
    1369 
    1370         /* Get node coordinates and dof list: */
    1371         GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
    1372         GetDofList(&doflist[0],&numberofdofspernode);
    1373 
    1374         /* Set pe_g to 0: */
    1375         for(i=0;i<numdof;i++) pe_g[i]=0.0;
    1376 
    1377         if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==SurfaceYAnalysisEnum())){
    1378                 for(i=0;i<numdof;i++) param[i]=s[i];
    1379         }
    1380         if ( (sub_analysis_type==BedXAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
    1381                 for(i=0;i<numdof;i++) param[i]=b[i];
    1382         }
    1383 
    1384         /* Get gaussian points and weights: */
    1385         GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/
    1386 
    1387 
    1388         /* Start  looping on the number of gaussian points: */
    1389         for (ig=0; ig<num_gauss; ig++){
    1390                 /*Pick up the gaussian point: */
    1391                 gauss_weight=*(gauss_weights+ig);
    1392                 gauss_l1l2l3[0]=*(first_gauss_area_coord+ig);
    1393                 gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
    1394                 gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
    1395 
    1396                 GetParameterDerivativeValue(&slope[0], &param[0],&xyz_list[0][0], gauss_l1l2l3);
    1397                
    1398                 /* Get Jacobian determinant: */
    1399                 GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
    1400                
    1401                  /*Get nodal functions: */
    1402                 GetNodalFunctions(l1l2l3, gauss_l1l2l3);
    1403 
    1404                 /*Build pe_g_gaussian vector: */
    1405                 if ( (sub_analysis_type==SurfaceXAnalysisEnum()) || (sub_analysis_type==BedXAnalysisEnum())){
    1406                         for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[0]*l1l2l3[i];
    1407                 }
    1408                 if ( (sub_analysis_type==SurfaceYAnalysisEnum()) || (sub_analysis_type==BedYAnalysisEnum())){
    1409                         for(i=0;i<numdof;i++) pe_g_gaussian[i]=Jdet*gauss_weight*slope[1]*l1l2l3[i];
    1410                 }
    1411 
    1412                 /*Add pe_g_gaussian vector to pe_g: */
    1413                 for( i=0; i<numdof; i++)pe_g[i]+=pe_g_gaussian[i];
    1414 
    1415         } //for (ig=0; ig<num_gauss; ig++)
    1416 
    1417         /*Add pe_g to global vector pg: */
    1418         VecSetValues(pg,numdof,doflist,(const double*)pe_g,ADD_VALUES);
    1419 
    1420         cleanup_and_return:
    1421         xfree((void**)&first_gauss_area_coord);
    1422         xfree((void**)&second_gauss_area_coord);
    1423         xfree((void**)&third_gauss_area_coord);
    1424         xfree((void**)&gauss_weights);
    1425 
    1426 }
    1427 
     1450}
     1451/*}}}*/
     1452/*FUNCTION UpdateFromInputs {{{1*/
    14281453#undef __FUNCT__
    14291454#define __FUNCT__ "Tria::UpdateFromInputs"
    … …  
    14631488
    14641489}
    1465                
     1490/*}}}*/
     1491/*FUNCTION GetDofList {{{1*/
    14661492void  Tria::GetDofList(int* doflist,int* pnumberofdofspernode){
    14671493
    … …  
    14811507
    14821508}
    1483 
     1509/*}}}*/
     1510/*FUNCTION GetDofList1 {{{1*/
    14841511void  Tria::GetDofList1(int* doflist){
    14851512
    … …  
    14901517
    14911518}
    1492 
    1493 
     1519/*}}}*/
     1520/*FUNCTION GetParameterValue {{{1*/
    14941521#undef __FUNCT__
    14951522#define __FUNCT__ "Tria::GetParameterValue"
    … …  
    15121539        *pp=p;
    15131540}
    1514 
    1515 
     1541/*}}}*/
     1542/*FUNCTION GetParameterDerivativeValue {{{1*/
    15161543#undef __FUNCT__
    15171544#define __FUNCT__ "Tria::GetParameterDerivativeValue"
    … …  
    15371564
    15381565}
    1539 
    1540 
     1566/*}}}*/
     1567/*FUNCTION GetStrainRate {{{1*/
    15411568#undef __FUNCT__
    15421569#define __FUNCT__ "Tria::GetStrainRate"
    … …  
    15581585
    15591586}
    1560 
     1587/*}}}*/
     1588/*FUNCTION GetJacobianDeterminant2d {{{1*/
    15611589#undef __FUNCT__
    15621590#define __FUNCT__ "Tria::GetJacobianDeterminant2d"
    … …  
    15841612       
    15851613}
    1586 
     1614/*}}}*/
     1615/*FUNCTION GetJacobianDeterminant3d {{{1*/
    15871616#undef __FUNCT__
    15881617#define __FUNCT__ "Tria::GetJacobianDeterminant3d"
    … …  
    16131642       
    16141643}
    1615 
     1644/*}}}*/
     1645/*FUNCTION GetB {{{1*/
    16161646#undef __FUNCT__
    16171647#define __FUNCT__ "Tria::GetB"
    … …  
    16561686        }
    16571687}
    1658 
     1688/*}}}*/
     1689/*FUNCTION GetBPrime {{{1*/
    16591690#undef __FUNCT__
    16601691#define __FUNCT__ "Tria::GetBPrime"
    … …  
    16941725        }
    16951726}
    1696 
     1727/*}}}*/
     1728/*FUNCTION GetL {{{1*/
    16971729#undef __FUNCT__
    16981730#define __FUNCT__ "Tria::GetL"
    … …  
    17441776        }
    17451777}
    1746 
     1778/*}}}*/
     1779/*FUNCTION GetB_prog {{{1*/
    17471780#undef __FUNCT__
    17481781#define __FUNCT__ "Tria::GetB_prog"
    … …  
    17821815        }
    17831816}
    1784 
     1817/*}}}*/
     1818/*FUNCTION GetBPrime_prog {{{1*/
    17851819#undef __FUNCT__
    17861820#define __FUNCT__ "Tria::GetBPrime_prog"
    … …  
    18151849        }
    18161850}
    1817 
    1818 
     1851/*}}}*/
     1852/*FUNCTION GetNodalFunctions {{{1*/
    18191853#undef __FUNCT__
    18201854#define __FUNCT__ "Tria::GetNodalFunctions"
    … …  
    18331867
    18341868}
    1835 
     1869/*}}}*/
     1870/*FUNCTION GetNodalFunctionsDerivativesBasic {{{1*/
    18361871#undef __FUNCT__
    18371872#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesBasic"
    … …  
    18671902
    18681903}
    1869 
     1904/*}}}*/
     1905/*FUNCTION GetNodalFunctionsDerivativesParams {{{1*/
    18701906#undef __FUNCT__
    18711907#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesParams"
    … …  
    18931929
    18941930}
    1895 
     1931/*}}}*/
     1932/*FUNCTION GetJacobianInvert {{{1*/
    18961933#undef __FUNCT__
    18971934#define __FUNCT__ "Tria::GetJacobianInvert"
    … …  
    19091946               
    19101947}
    1911 
     1948/*}}}*/
     1949/*FUNCTION GetJacobian {{{1*/
    19121950#undef __FUNCT__
    19131951#define __FUNCT__ "Tria::GetJacobian"
    … …  
    19772015
    19782016}
    1979 
    1980 
     2017/*}}}*/
     2018/*FUNCTION Du {{{1*/
    19812019#undef __FUNCT__
    19822020#define __FUNCT__ "Tria::Du"
    … …  
    21892227
    21902228}
     2229/*}}}*/
     2230/*FUNCTION Gradj {{{1*/
    21912231#undef __FUNCT__
    21922232#define __FUNCT__ "Tria::Gradj"
    … …  
    22042244        else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
    22052245}
    2206 
     2246/*}}}*/
     2247/*FUNCTION GradjDrag {{{1*/
    22072248#undef __FUNCT__
    22082249#define __FUNCT__ "Tria::GradjDrag"
    … …  
    24232464
    24242465}
    2425 
     2466/*}}}*/
     2467/*FUNCTION GradjDragStokes {{{1*/
    24262468#undef __FUNCT__
    24272469#define __FUNCT__ "Tria::GradjDragStokes"
    … …  
    26582700
    26592701}
    2660 
     2702/*}}}*/
     2703/*FUNCTION SurfaceNormal{{{1*/
    26612704#undef __FUNCT__
    26622705#define __FUNCT__ "Tria::SurfaceNormal"
    … …  
    26862729
    26872730}
    2688 
     2731/*}}}*/
     2732/*FUNCTION GradjB{{{1*/
    26892733#undef __FUNCT__
    26902734#define __FUNCT__ "Tria::GradjB"
    … …  
    28612905        xfree((void**)&gauss_weights);
    28622906}
    2863 
    2864 
     2907/*}}}*/
     2908/*FUNCTION Misfit {{{1*/
    28652909#undef __FUNCT__
    28662910#define __FUNCT__ "Tria::Misfit"
    … …  
    30653109        return Jelem;
    30663110}
    3067 
     3111/*}}}*/
     3112/*FUNCTION NodeConfiguration {{{1*/
    30683113#undef __FUNCT__
    30693114#define __FUNCT__ "Tria::NodeConfiguration"
    … …  
    30783123
    30793124}
     3125/*}}}*/
     3126/*FUNCTION MaticeConfiguration {{{1*/
    30803127#undef __FUNCT__
    30813128#define __FUNCT__ "Tria::MaticeConfiguration"
    … …  
    30843131        matice_offset=tria_matice_offset;
    30853132}
    3086 
    3087 #undef __FUNCT__
    3088 #define __FUNCT__ "Tria::MaticeConfiguration"
     3133/*}}}*/
     3134/*FUNCTION MatparConfiguration {{{1*/
     3135#undef __FUNCT__
     3136#define __FUNCT__ "Tria::MatparConfiguration"
    30893137void  Tria::MatparConfiguration(Matpar* tria_matpar,int tria_matpar_offset){
    30903138
    … …  
    30933141
    30943142}
    3095 
     3143/*}}}*/
     3144/*FUNCTION NumparConfiguration {{{1*/
    30963145#undef __FUNCT__
    30973146#define __FUNCT__ "Tria::NumparConfiguration"
    … …  
    31023151
    31033152}
    3104 
     3153/*}}}*/
     3154/*FUNCTION CreateKMatrixDiagnosticSurfaceVert {{{1*/
    31053155#undef __FUNCT__
    31063156#define __FUNCT__ "Tria::CreateKMatrixDiagnosticSurfaceVert"
    … …  
    32293279        xfree((void**)&gauss_weights);
    32303280}
    3231                
    3232 
     3281/*}}}*/
     3282/*FUNCTION CreatePVectorDiagnosticBaseVert {{{1*/
    32333283#undef __FUNCT__
    32343284#define __FUNCT__ "Tria::CreatePVectorDiagnosticBaseVert"
    … …  
    33463396
    33473397}
    3348 
     3398/*}}}*/
     3399/*FUNCTION ComputePressure {{{1*/
    33493400#undef __FUNCT__
    33503401#define __FUNCT__ "Tria::ComputePressure"
    … …  
    33713422
    33723423}
    3373 
     3424/*}}}*/
     3425/*FUNCTION CreateKMatrixThermal {{{1*/
    33743426#undef __FUNCT__
    33753427#define __FUNCT__ "Tria::CreateKMatrixThermal"
    … …  
    34723524
    34733525}
    3474 
    3475 
     3526/*}}}*/
     3527/*FUNCTION CreateKMatrixMelting {{{1*/
    34763528#undef __FUNCT__
    34773529#define __FUNCT__ "Tria::CreateKMatrixMelting"
    … …  
    35583610
    35593611}
    3560 
    3561 
     3612/*}}}*/
     3613/*FUNCTION CreatePVectorThermalShelf {{{1*/
    35623614#undef __FUNCT__
    35633615#define __FUNCT__ "Tria::CreatePVectorThermalShelf"
    … …  
    36713723
    36723724}
     3725/*}}}*/
     3726/*FUNCTION CreatePVectorThermalSheet {{{1*/
    36733727#undef __FUNCT__
    36743728#define __FUNCT__ "Tria::CreatePVectorThermalSheet"
    … …  
    38143868
    38153869}
    3816 
     3870/*}}}*/
     3871/*FUNCTION MassFlux {{{1*/
    38173872#undef __FUNCT__
    38183873#define __FUNCT__ "Tria::MassFlux"
    … …  
    38853940        return mass_flux;
    38863941}
    3887 
     3942/*}}}*/
     3943/*FUNCTION GetArea {{{1*/
    38883944#undef __FUNCT__
    38893945#define __FUNCT__ "Tria::GetArea"
    … …  
    39033959        return x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1;
    39043960}
    3905 
     3961/*}}}*/
     3962/*FUNCTION GetAreaCoordinate {{{1*/
    39063963#undef __FUNCT__
    39073964#define __FUNCT__ "Tria::GetAreaCoordinate"
    … …  
    39363993        else throw ErrorException(__FUNCT__,exprintf("%s%i%s\n"," error message: area coordinate ",which_one," done not exist!"));
    39373994}
    3938 
    3939 
     3995/*}}}*/
Note: See TracChangeset for help on using the changeset viewer.