Ignore:
Timestamp:
02/13/12 15:50:19 (13 years ago)
Author:
cborstad
Message:

merged src changes 11330:11410 from trunk-jpl

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl-damage/src/c/objects/Elements/Penta.cpp

    r11409 r11417  
    769769}
    770770/*}}}*/
     771<<<<<<< .working
     772=======
     773/*FUNCTION Penta::CreateJacobianMatrix{{{1*/
     774void  Penta::CreateJacobianMatrix(Mat Jff){
     775
     776        /*retrieve parameters: */
     777        ElementMatrix* Ke=NULL;
     778        int analysis_type;
     779        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     780
     781        /*Checks in debugging {{{2*/
     782        _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
     783        /*}}}*/
     784
     785        /*Skip if water element*/
     786        if(IsOnWater()) return;
     787
     788        /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
     789        switch(analysis_type){
     790#ifdef _HAVE_DIAGNOSTIC_
     791                case DiagnosticHorizAnalysisEnum:
     792                        Ke=CreateJacobianDiagnosticHoriz();
     793                        break;
     794#endif
     795                default:
     796                        _error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
     797        }
     798
     799        /*Add to global matrix*/
     800        if(Ke){
     801                Ke->AddToGlobal(Jff);
     802                delete Ke;
     803        }
     804}
     805/*}}}*/
     806>>>>>>> .merge-right.r11410
    771807/*FUNCTION Penta::DeepEcho{{{1*/
    772808void Penta::DeepEcho(void){
    … …  
    10991135/*}}}*/
    11001136/*FUNCTION Penta::GetStabilizationParameter {{{1*/
    1101 double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity){
     1137double Penta::GetStabilizationParameter(double u, double v, double w, double diameter, double kappa){
    11021138        /*Compute stabilization parameter*/
     1139        /*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
     1140        /*kappa=enthalpydiffusionparameter for enthalpy model*/
    11031141
    11041142        double normu;
    … …  
    11061144
    11071145        normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
    1108         if(normu*diameter/(3*2*thermalconductivity/(rho_ice*heatcapacity))<1){
    1109                 tau_parameter=pow(diameter,2)/(3*2*2*thermalconductivity/(rho_ice*heatcapacity));
     1146        if(normu*diameter/(3*2*kappa)<1){
     1147                tau_parameter=pow(diameter,2)/(3*2*2*kappa);
    11101148        }
    11111149        else tau_parameter=diameter/(2*normu);
    … …  
    32833321                GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss);
    32843322
    3285                 vx_input->GetInputValue(&u, gauss);
    3286                 vy_input->GetInputValue(&v, gauss);
    3287                 vz_input->GetInputValue(&w, gauss);
    3288                 vxm_input->GetInputValue(&um,gauss);
    3289                 vym_input->GetInputValue(&vm,gauss);
    3290                 vzm_input->GetInputValue(&wm,gauss);
    3291                 vx=u-um; vy=v-vm; vz=w-wm;
     3323                vx_input->GetInputValue(&u, gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
     3324                vy_input->GetInputValue(&v, gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
     3325                vz_input->GetInputValue(&w, gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    32923326
    32933327                D_scalar_advec=gauss->weight*Jdet;
    … …  
    33213355                        vel=sqrt(pow(vx,2.)+pow(vy,2.)+pow(vz,2.))+1.e-14;
    33223356                        h=sqrt( pow(hx*vx/vel,2.) + pow(hy*vy/vel,2.) + pow(hz*vz/vel,2.));
    3323                         K[0][0]=h/(2*vel)*fabs(vx*vx);  K[0][1]=h/(2*vel)*fabs(vx*vy); K[0][2]=h/(2*vel)*fabs(vx*vz);
    3324                         K[1][0]=h/(2*vel)*fabs(vy*vx);  K[1][1]=h/(2*vel)*fabs(vy*vy); K[1][2]=h/(2*vel)*fabs(vy*vz);
    3325                         K[2][0]=h/(2*vel)*fabs(vz*vx);  K[2][1]=h/(2*vel)*fabs(vz*vy); K[2][2]=h/(2*vel)*fabs(vz*vz);
     3357                        K[0][0]=h/(2*vel)*vx*vx;  K[0][1]=h/(2*vel)*vx*vy; K[0][2]=h/(2*vel)*vx*vz;
     3358                        K[1][0]=h/(2*vel)*vy*vx;  K[1][1]=h/(2*vel)*vy*vy; K[1][2]=h/(2*vel)*vy*vz;
     3359                        K[2][0]=h/(2*vel)*vz*vx;  K[2][1]=h/(2*vel)*vz*vy; K[2][2]=h/(2*vel)*vz*vz;
    33263360                        D_scalar_stab=gauss->weight*Jdet;
    33273361                        if(dt) D_scalar_stab=D_scalar_stab*dt;
    … …  
    33373371                else if(stabilization==2){
    33383372                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3339                         tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
     3373                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
    33403374
    33413375                        for(i=0;i<numdof;i++){
    … …  
    34513485        double     Jdet,u,v,w,um,vm,wm,vel;
    34523486        double     h,hx,hy,hz,vx,vy,vz;
    3453         double     gravity,rho_ice,rho_water;
     3487        double     gravity,rho_ice,rho_water,kappa;
    34543488        double     heatcapacity,thermalconductivity,dt;
    34553489        double     tau_parameter,diameter;
    … …  
    34773511        heatcapacity=matpar->GetHeatCapacity();
    34783512        thermalconductivity=matpar->GetThermalConductivity();
     3513        kappa=thermalconductivity/(rho_ice*heatcapacity);
    34793514        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    34803515        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    … …  
    34993534                GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss);
    35003535
    3501                 D_scalar_conduct=gauss->weight*Jdet*(thermalconductivity/(rho_ice*heatcapacity));
     3536                D_scalar_conduct=gauss->weight*Jdet*kappa;
    35023537                if(dt) D_scalar_conduct=D_scalar_conduct*dt;
    35033538
    … …  
    35683603                else if(stabilization==2){
    35693604                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3570                         tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
     3605                        tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,kappa);
    35713606
    35723607                        for(i=0;i<numdof;i++){
    … …  
    36733708        double Jdet,phi,dt;
    36743709        double rho_ice,heatcapacity;
    3675         double thermalconductivity;
    3676         double viscosity,enthalpy;
     3710        double thermalconductivity,kappa;
     3711        double viscosity,pressure;
     3712        double enthalpy,enthalpypicard;
    36773713        double tau_parameter,diameter;
    36783714        double u,v,w;
    … …  
    36953731        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    36963732        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    3697         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    3698         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    3699         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    3700         Input* enthalpy_input=NULL;
    3701         if (dt) enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(inputs);
    3702         if (stabilization==2) diameter=MinEdgeLength(xyz_list);
     3733        Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
     3734        Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
     3735        Input* vz_input=inputs->GetInput(VzEnum);                                  _assert_(vz_input);
     3736        Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
     3737        Input* enthalpy_input=NULL;
     3738        Input* enthalpypicard_input=NULL;
     3739        if(dt){
     3740                enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
     3741        }
     3742        if (stabilization==2){
     3743                diameter=MinEdgeLength(xyz_list);
     3744                enthalpypicard_input=inputs->GetInput(EnthalpyPicardEnum); _assert_(enthalpypicard_input);
     3745        }
    37033746
    37043747        /* Start  looping on the number of gaussian points: */
    … …  
    37153758                GetPhi(&phi, &epsilon[0], viscosity);
    37163759
    3717                 scalar_def=phi/(rho_ice)*Jdet*gauss->weight;
     3760                scalar_def=phi/rho_ice*Jdet*gauss->weight;
    37183761                if(dt) scalar_def=scalar_def*dt;
    37193762
    … …  
    37333776                        vy_input->GetInputValue(&v, gauss);
    37343777                        vz_input->GetInputValue(&w, gauss);
    3735 
    3736                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
     3778                        pressure_input->GetInputValue(&pressure, gauss);
     3779                        enthalpypicard_input->GetInputValue(&enthalpypicard, gauss);
     3780                        kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
     3781                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
    37373782
    37383783                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    … …  
    38193864        double     rho_ice,heatcapacity,geothermalflux_value;
    38203865        double     basalfriction,alpha2,vx,vy;
     3866        double     scalar,enthalpy,enthalpyup;
     3867        double     pressure,pressureup;
     3868        double     basis[NUMVERTICES];
     3869        Friction*  friction=NULL;
     3870        GaussPenta* gauss=NULL;
     3871        GaussPenta* gaussup=NULL;
     3872
     3873        /* Geothermal flux on ice sheet base and basal friction */
     3874        if (!IsOnBed() || IsFloating()) return NULL;
     3875
     3876        /*Initialize Element vector*/
     3877        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3878
     3879        /*Retrieve all inputs and parameters*/
     3880        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3881        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     3882        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
     3883        rho_ice=matpar->GetRhoIce();
     3884        heatcapacity=matpar->GetHeatCapacity();
     3885        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     3886        Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
     3887        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
     3888        Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
     3889        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
     3890        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
     3891        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
     3892
     3893        /*Build frictoin element, needed later: */
     3894        friction=new Friction("3d",inputs,matpar,analysis_type);
     3895
     3896        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     3897        gauss=new GaussPenta(0,1,2,2);
     3898        gaussup=new GaussPenta(3,4,5,2);
     3899        for(ig=gauss->begin();ig<gauss->end();ig++){
     3900
     3901                gauss->GaussPoint(ig);
     3902                gaussup->GaussPoint(ig);
     3903
     3904                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     3905                GetNodalFunctionsP1(&basis[0], gauss);
     3906
     3907                enthalpy_input->GetInputValue(&enthalpy,gauss);
     3908                pressure_input->GetInputValue(&pressure,gauss);
     3909//              if(enthalpy>matpar->PureIceEnthalpy(pressure)){
     3910//                      enthalpy_input->GetInputValue(&enthalpyup,gaussup);
     3911//                      pressure_input->GetInputValue(&pressureup,gaussup);
     3912//                      if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
     3913//                              //do nothing, don't add heatflux
     3914//                      }
     3915//                      else{
     3916//                              //need to change spcenthalpy according to Aschwanden
     3917//                              //NEED TO UPDATE
     3918//                      }
     3919//              }
     3920//              else{
     3921                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     3922                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     3923                        vx_input->GetInputValue(&vx,gauss);
     3924                        vy_input->GetInputValue(&vy,gauss);
     3925                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     3926
     3927                        scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
     3928                        if(dt) scalar=dt*scalar;
     3929
     3930                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     3931//              }
     3932        }
     3933
     3934        /*Clean up and return*/
     3935        delete gauss;
     3936        delete gaussup;
     3937        delete friction;
     3938        return pe;
     3939}
     3940/*}}}*/
     3941/*FUNCTION Penta::CreatePVectorMelting {{{1*/
     3942ElementVector* Penta::CreatePVectorMelting(void){
     3943        return NULL;
     3944}
     3945/*}}}*/
     3946/*FUNCTION Penta::CreatePVectorThermal {{{1*/
     3947ElementVector* Penta::CreatePVectorThermal(void){
     3948
     3949        /*compute all load vectors for this element*/
     3950        ElementVector* pe1=CreatePVectorThermalVolume();
     3951        ElementVector* pe2=CreatePVectorThermalSheet();
     3952        ElementVector* pe3=CreatePVectorThermalShelf();
     3953        ElementVector* pe =new ElementVector(pe1,pe2,pe3);
     3954
     3955        /*clean-up and return*/
     3956        delete pe1;
     3957        delete pe2;
     3958        delete pe3;
     3959        return pe;
     3960}
     3961/*}}}*/
     3962/*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/
     3963ElementVector* Penta::CreatePVectorThermalVolume(void){
     3964
     3965        /*Constants*/
     3966        const int  numdof=NUMVERTICES*NDOF1;
     3967
     3968        /*Intermediaries*/
     3969        int    i,j,ig,found=0;
     3970        int    friction_type,stabilization;
     3971        double Jdet,phi,dt;
     3972        double rho_ice,heatcapacity;
     3973        double thermalconductivity,kappa;
     3974        double viscosity,temperature;
     3975        double tau_parameter,diameter;
     3976        double u,v,w;
     3977        double scalar_def,scalar_transient;
     3978        double temperature_list[NUMVERTICES];
     3979        double xyz_list[NUMVERTICES][3];
     3980        double L[numdof];
     3981        double dbasis[3][6];
     3982        double epsilon[6];
     3983        GaussPenta *gauss=NULL;
     3984
     3985        /*Initialize Element vector*/
     3986        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     3987
     3988        /*Retrieve all inputs and parameters*/
     3989        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     3990        rho_ice=matpar->GetRhoIce();
     3991        heatcapacity=matpar->GetHeatCapacity();
     3992        thermalconductivity=matpar->GetThermalConductivity();
     3993        kappa=thermalconductivity/(rho_ice*heatcapacity);
     3994        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     3995        this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
     3996        Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
     3997        Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
     3998        Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
     3999        Input* temperature_input=NULL;
     4000        if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
     4001        if (stabilization==2) diameter=MinEdgeLength(xyz_list);
     4002
     4003        /* Start  looping on the number of gaussian points: */
     4004        gauss=new GaussPenta(2,3);
     4005        for (ig=gauss->begin();ig<gauss->end();ig++){
     4006
     4007                gauss->GaussPoint(ig);
     4008
     4009                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     4010                GetNodalFunctionsP1(&L[0], gauss);
     4011
     4012                this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     4013                matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
     4014                GetPhi(&phi, &epsilon[0], viscosity);
     4015
     4016                scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
     4017                if(dt) scalar_def=scalar_def*dt;
     4018
     4019                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
     4020
     4021                /* Build transient now */
     4022                if(dt){
     4023                        temperature_input->GetInputValue(&temperature, gauss);
     4024                        scalar_transient=temperature*Jdet*gauss->weight;
     4025                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
     4026                }
     4027
     4028                if(stabilization==2){
     4029                        GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
     4030
     4031                        vx_input->GetInputValue(&u, gauss);
     4032                        vy_input->GetInputValue(&v, gauss);
     4033                        vz_input->GetInputValue(&w, gauss);
     4034
     4035                        tau_parameter=GetStabilizationParameter(u,v,w,diameter,kappa);
     4036
     4037                        for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
     4038                        if(dt){
     4039                                for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
     4040                        }
     4041                }
     4042        }
     4043
     4044        /*Clean up and return*/
     4045        delete gauss;
     4046        return pe;
     4047}
     4048/*}}}*/
     4049/*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/
     4050ElementVector* Penta::CreatePVectorThermalShelf(void){
     4051
     4052        /*Constants*/
     4053        const int  numdof=NUMVERTICES*NDOF1;
     4054
     4055        /*Intermediaries */
     4056        int        i,j,ig;
     4057        double     Jdet2d;
     4058        double     mixed_layer_capacity,thermal_exchange_velocity;
     4059        double     rho_ice,rho_water,pressure,dt,scalar_ocean;
     4060        double     heatcapacity,t_pmp;
     4061        double     xyz_list[NUMVERTICES][3];
     4062        double     xyz_list_tria[NUMVERTICES2D][3];
     4063        double     basis[NUMVERTICES];
     4064        GaussPenta* gauss=NULL;
     4065
     4066        /* Ice/ocean heat exchange flux on ice shelf base */
     4067        if (!IsOnBed() || !IsFloating()) return NULL;
     4068
     4069        /*Initialize Element vector*/
     4070        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     4071
     4072        /*Retrieve all inputs and parameters*/
     4073        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     4074        for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
     4075        mixed_layer_capacity=matpar->GetMixedLayerCapacity();
     4076        thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
     4077        rho_water=matpar->GetRhoWater();
     4078        rho_ice=matpar->GetRhoIce();
     4079        heatcapacity=matpar->GetHeatCapacity();
     4080        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4081        Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
     4082
     4083        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
     4084        gauss=new GaussPenta(0,1,2,2);
     4085        for(ig=gauss->begin();ig<gauss->end();ig++){
     4086
     4087                gauss->GaussPoint(ig);
     4088
     4089                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
     4090                GetNodalFunctionsP1(&basis[0], gauss);
     4091
     4092                pressure_input->GetInputValue(&pressure,gauss);
     4093                t_pmp=matpar->TMeltingPoint(pressure);
     4094
     4095                scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
     4096                if(dt) scalar_ocean=dt*scalar_ocean;
     4097
     4098                for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
     4099        }
     4100
     4101        /*Clean up and return*/
     4102        delete gauss;
     4103        return pe;
     4104}
     4105/*}}}*/
     4106/*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/
     4107ElementVector* Penta::CreatePVectorThermalSheet(void){
     4108
     4109        /*Constants*/
     4110        const int  numdof=NUMVERTICES*NDOF1;
     4111
     4112        /*Intermediaries */
     4113        int        i,j,ig;
     4114        int        analysis_type;
     4115        double     xyz_list[NUMVERTICES][3];
     4116        double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
     4117        double     Jdet2d,dt;
     4118        double     rho_ice,heatcapacity,geothermalflux_value;
     4119        double     basalfriction,alpha2,vx,vy;
     4120        double     basis[NUMVERTICES];
    38214121        double     scalar;
    3822         double     basis[NUMVERTICES];
    38234122        Friction*  friction=NULL;
    38244123        GaussPenta* gauss=NULL;
    … …  
    38544153                GetNodalFunctionsP1(&basis[0], gauss);
    38554154
    3856                 geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    3857                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    3858                 vx_input->GetInputValue(&vx,gauss);
    3859                 vy_input->GetInputValue(&vy,gauss);
    3860                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    3861                
    3862                 scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice);
    3863                 if(dt) scalar=dt*scalar;
    3864 
    3865                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    3866         }
    3867 
    3868         /*Clean up and return*/
    3869         delete gauss;
    3870         delete friction;
    3871         return pe;
    3872 }
    3873 /*}}}*/
    3874 /*FUNCTION Penta::CreatePVectorMelting {{{1*/
    3875 ElementVector* Penta::CreatePVectorMelting(void){
    3876         return NULL;
    3877 }
    3878 /*}}}*/
    3879 /*FUNCTION Penta::CreatePVectorThermal {{{1*/
    3880 ElementVector* Penta::CreatePVectorThermal(void){
    3881 
    3882         /*compute all load vectors for this element*/
    3883         ElementVector* pe1=CreatePVectorThermalVolume();
    3884         ElementVector* pe2=CreatePVectorThermalSheet();
    3885         ElementVector* pe3=CreatePVectorThermalShelf();
    3886         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    3887 
    3888         /*clean-up and return*/
    3889         delete pe1;
    3890         delete pe2;
    3891         delete pe3;
    3892         return pe;
    3893 }
    3894 /*}}}*/
    3895 /*FUNCTION Penta::CreatePVectorThermalVolume {{{1*/
    3896 ElementVector* Penta::CreatePVectorThermalVolume(void){
    3897 
    3898         /*Constants*/
    3899         const int  numdof=NUMVERTICES*NDOF1;
    3900 
    3901         /*Intermediaries*/
    3902         int    i,j,ig,found=0;
    3903         int    friction_type,stabilization;
    3904         double Jdet,phi,dt;
    3905         double rho_ice,heatcapacity;
    3906         double thermalconductivity;
    3907         double viscosity,temperature;
    3908         double tau_parameter,diameter;
    3909         double u,v,w;
    3910         double scalar_def,scalar_transient;
    3911         double temperature_list[NUMVERTICES];
    3912         double xyz_list[NUMVERTICES][3];
    3913         double L[numdof];
    3914         double dbasis[3][6];
    3915         double epsilon[6];
    3916         GaussPenta *gauss=NULL;
    3917 
    3918         /*Initialize Element vector*/
    3919         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    3920 
    3921         /*Retrieve all inputs and parameters*/
    3922         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    3923         rho_ice=matpar->GetRhoIce();
    3924         heatcapacity=matpar->GetHeatCapacity();
    3925         thermalconductivity=matpar->GetThermalConductivity();
    3926         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    3927         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    3928         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    3929         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    3930         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    3931         Input* temperature_input=NULL;
    3932         if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
    3933         if (stabilization==2) diameter=MinEdgeLength(xyz_list);
    3934 
    3935         /* Start  looping on the number of gaussian points: */
    3936         gauss=new GaussPenta(2,3);
    3937         for (ig=gauss->begin();ig<gauss->end();ig++){
    3938 
    3939                 gauss->GaussPoint(ig);
    3940 
    3941                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3942                 GetNodalFunctionsP1(&L[0], gauss);
    3943 
    3944                 this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    3945                 matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
    3946                 GetPhi(&phi, &epsilon[0], viscosity);
    3947 
    3948                 scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
    3949                 if(dt) scalar_def=scalar_def*dt;
    3950 
    3951                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
    3952 
    3953                 /* Build transient now */
    3954                 if(dt){
    3955                         temperature_input->GetInputValue(&temperature, gauss);
    3956                         scalar_transient=temperature*Jdet*gauss->weight;
    3957                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
    3958                 }
    3959 
    3960                 if(stabilization==2){
    3961                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    3962 
    3963                         vx_input->GetInputValue(&u, gauss);
    3964                         vy_input->GetInputValue(&v, gauss);
    3965                         vz_input->GetInputValue(&w, gauss);
    3966 
    3967                         tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
    3968 
    3969                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    3970                         if(dt){
    3971                                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
    3972                         }
    3973                 }
    3974         }
    3975 
    3976         /*Clean up and return*/
    3977         delete gauss;
    3978         return pe;
    3979 }
    3980 /*}}}*/
    3981 /*FUNCTION Penta::CreatePVectorThermalShelf {{{1*/
    3982 ElementVector* Penta::CreatePVectorThermalShelf(void){
    3983 
    3984         /*Constants*/
    3985         const int  numdof=NUMVERTICES*NDOF1;
    3986 
    3987         /*Intermediaries */
    3988         int        i,j,ig;
    3989         double     Jdet2d;
    3990         double     mixed_layer_capacity,thermal_exchange_velocity;
    3991         double     rho_ice,rho_water,pressure,dt,scalar_ocean;
    3992         double     heatcapacity,t_pmp;
    3993         double     xyz_list[NUMVERTICES][3];
    3994         double     xyz_list_tria[NUMVERTICES2D][3];
    3995         double     basis[NUMVERTICES];
    3996         GaussPenta* gauss=NULL;
    3997 
    3998         /* Ice/ocean heat exchange flux on ice shelf base */
    3999         if (!IsOnBed() || !IsFloating()) return NULL;
    4000 
    4001         /*Initialize Element vector*/
    4002         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4003 
    4004         /*Retrieve all inputs and parameters*/
    4005         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4006         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4007         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    4008         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    4009         rho_water=matpar->GetRhoWater();
    4010         rho_ice=matpar->GetRhoIce();
    4011         heatcapacity=matpar->GetHeatCapacity();
    4012         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4013         Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    4014 
    4015         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    4016         gauss=new GaussPenta(0,1,2,2);
    4017         for(ig=gauss->begin();ig<gauss->end();ig++){
    4018 
    4019                 gauss->GaussPoint(ig);
    4020 
    4021                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    4022                 GetNodalFunctionsP1(&basis[0], gauss);
    4023 
    4024                 pressure_input->GetInputValue(&pressure,gauss);
    4025                 t_pmp=matpar->TMeltingPoint(pressure);
    4026 
    4027                 scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
    4028                 if(dt) scalar_ocean=dt*scalar_ocean;
    4029 
    4030                 for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
    4031         }
    4032 
    4033         /*Clean up and return*/
    4034         delete gauss;
    4035         return pe;
    4036 }
    4037 /*}}}*/
    4038 /*FUNCTION Penta::CreatePVectorThermalSheet {{{1*/
    4039 ElementVector* Penta::CreatePVectorThermalSheet(void){
    4040 
    4041         /*Constants*/
    4042         const int  numdof=NUMVERTICES*NDOF1;
    4043 
    4044         /*Intermediaries */
    4045         int        i,j,ig;
    4046         int        analysis_type;
    4047         double     xyz_list[NUMVERTICES][3];
    4048         double     xyz_list_tria[NUMVERTICES2D][3]={0.0};
    4049         double     Jdet2d,dt;
    4050         double     rho_ice,heatcapacity,geothermalflux_value;
    4051         double     basalfriction,alpha2,vx,vy;
    4052         double     scalar;
    4053         double     basis[NUMVERTICES];
    4054         Friction*  friction=NULL;
    4055         GaussPenta* gauss=NULL;
    4056 
    4057         /* Geothermal flux on ice sheet base and basal friction */
    4058         if (!IsOnBed() || IsFloating()) return NULL;
    4059 
    4060         /*Initialize Element vector*/
    4061         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4062 
    4063         /*Retrieve all inputs and parameters*/
    4064         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    4065         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4066         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4067         rho_ice=matpar->GetRhoIce();
    4068         heatcapacity=matpar->GetHeatCapacity();
    4069         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4070         Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
    4071         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
    4072         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
    4073         Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    4074 
    4075         /*Build frictoin element, needed later: */
    4076         friction=new Friction("3d",inputs,matpar,analysis_type);
    4077 
    4078         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    4079         gauss=new GaussPenta(0,1,2,2);
    4080         for(ig=gauss->begin();ig<gauss->end();ig++){
    4081 
    4082                 gauss->GaussPoint(ig);
    4083 
    4084                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    4085                 GetNodalFunctionsP1(&basis[0], gauss);
    4086 
    4087                 geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4088                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    4089                 vx_input->GetInputValue(&vx,gauss);
    4090                 vy_input->GetInputValue(&vy,gauss);
    4091                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    4092                
    4093                 scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    4094                 if(dt) scalar=dt*scalar;
    4095 
    4096                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     4155                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     4156                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4157                        vx_input->GetInputValue(&vx,gauss);
     4158                        vy_input->GetInputValue(&vy,gauss);
     4159                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     4160
     4161                        scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
     4162                        if(dt) scalar=dt*scalar;
     4163
     4164                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    40974165        }
    40984166
    … …  
    42714339        if(converged){
    42724340                /*Convert enthalpy into temperature and water fraction*/
     4341<<<<<<< .working
    42734342                for(i=0;i<numdof;i++) matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
     4343=======
     4344                for(i=0;i<numdof;i++){
     4345                        matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
     4346                        if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector");
     4347                        //if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector");
     4348                }
     4349>>>>>>> .merge-right.r11410
    42744350                       
    42754351                this->inputs->AddInput(new PentaP1Input(EnthalpyEnum,values));
    … …  
    73197395}
    73207396/*}}}*/
     7397/*FUNCTION Penta::CreateJacobianDiagnosticHoriz{{{1*/
     7398ElementMatrix* Penta::CreateJacobianDiagnosticHoriz(void){
     7399
     7400        int approximation;
     7401        inputs->GetInputValue(&approximation,ApproximationEnum);
     7402
     7403        switch(approximation){
     7404                case MacAyealApproximationEnum:
     7405                        return CreateJacobianDiagnosticMacayeal2d();
     7406                case PattynApproximationEnum:
     7407                        return CreateJacobianDiagnosticPattyn();
     7408                case StokesApproximationEnum:
     7409                        return CreateJacobianDiagnosticStokes();
     7410                case NoneApproximationEnum:
     7411                        return NULL;
     7412                default:
     7413                        _error_("Approximation %s not supported yet",EnumToStringx(approximation));
     7414        }
     7415}
     7416/*}}}*/
     7417/*FUNCTION Penta::CreateJacobianDiagnosticMacayeal2d{{{1*/
     7418ElementMatrix* Penta::CreateJacobianDiagnosticMacayeal2d(void){
     7419
     7420        /*Figure out if this penta is collapsed. If so, then bailout, except if it is at the
     7421          bedrock, in which case we spawn a tria element using the 3 first nodes, and use it to build
     7422          the stiffness matrix. */
     7423        if (!IsOnBed()) return NULL;
     7424
     7425        /*Depth Averaging B*/
     7426        this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum,MaterialsEnum);
     7427
     7428        /*Call Tria function*/
     7429        Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
     7430        ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal();
     7431        delete tria->matice; delete tria;
     7432
     7433        /*Delete B averaged*/
     7434        this->matice->inputs->DeleteInput(MaterialsRheologyBbarEnum);
     7435
     7436        /*clean up and return*/
     7437        return Ke;
     7438}
     7439/*}}}*/
     7440/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{1*/
     7441ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
     7442
     7443        /*Constants*/
     7444        const int    numdof=NDOF2*NUMVERTICES;
     7445
     7446        /*Intermediaries */
     7447        int        i,j,ig;
     7448        double     xyz_list[NUMVERTICES][3];
     7449        double     Jdet;
     7450        double     eps1dotdphii,eps1dotdphij;
     7451        double     eps2dotdphii,eps2dotdphij;
     7452        double     mu_prime;
     7453        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     7454        double     eps1[3],eps2[3];
     7455        double     phi[NUMVERTICES];
     7456        double     dphi[3][NUMVERTICES];
     7457        GaussPenta *gauss=NULL;
     7458
     7459        /*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
     7460        ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
     7461
     7462        /*Retrieve all inputs and parameters*/
     7463        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     7464        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     7465        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     7466
     7467        /* Start  looping on the number of gaussian points: */
     7468        gauss=new GaussPenta(5,5);
     7469        for (ig=gauss->begin();ig<gauss->end();ig++){
     7470
     7471                gauss->GaussPoint(ig);
     7472
     7473                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7474                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
     7475
     7476                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     7477                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     7478                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
     7479                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
     7480                eps1[2]=epsilon[3];                eps2[2]=epsilon[4];
     7481
     7482                for(i=0;i<6;i++){
     7483                        for(j=0;j<6;j++){
     7484                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
     7485                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
     7486                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
     7487                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
     7488
     7489                                Ke->values[12*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     7490                                Ke->values[12*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     7491                                Ke->values[12*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     7492                                Ke->values[12*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     7493                        }
     7494                }
     7495        }
     7496
     7497        /*Transform Coordinate System*/
     7498        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);
     7499
     7500        /*Clean up and return*/
     7501        delete gauss;
     7502        return Ke;
     7503}
     7504/*}}}*/
     7505/*FUNCTION Penta::CreateJacobianDiagnosticStokes{{{1*/
     7506ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
     7507
     7508        /*Constants*/
     7509        const int    numdof=NDOF4*NUMVERTICES;
     7510
     7511        /*Intermediaries */
     7512        int        i,j,ig;
     7513        double     xyz_list[NUMVERTICES][3];
     7514        double     Jdet;
     7515        double     eps1dotdphii,eps1dotdphij;
     7516        double     eps2dotdphii,eps2dotdphij;
     7517        double     eps3dotdphii,eps3dotdphij;
     7518        double     mu_prime;
     7519        double     epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
     7520        double     eps1[3],eps2[3],eps3[3];
     7521        double     phi[NUMVERTICES];
     7522        double     dphi[3][NUMVERTICES];
     7523        GaussPenta *gauss=NULL;
     7524
     7525        /*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
     7526        ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
     7527
     7528        /*Retrieve all inputs and parameters*/
     7529        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     7530        Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
     7531        Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
     7532        Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
     7533
     7534        /* Start  looping on the number of gaussian points: */
     7535        gauss=new GaussPenta(5,5);
     7536        for (ig=gauss->begin();ig<gauss->end();ig++){
     7537
     7538                gauss->GaussPoint(ig);
     7539
     7540                GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
     7541                GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
     7542
     7543                this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
     7544                matice->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     7545                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
     7546                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
     7547                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
     7548
     7549                for(i=0;i<6;i++){
     7550                        for(j=0;j<6;j++){
     7551                                eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i]+eps1[2]*dphi[2][i];
     7552                                eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j]+eps1[2]*dphi[2][j];
     7553                                eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i]+eps2[2]*dphi[2][i];
     7554                                eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j]+eps2[2]*dphi[2][j];
     7555                                eps3dotdphii=eps3[0]*dphi[0][i]+eps3[1]*dphi[1][i]+eps3[2]*dphi[2][i];
     7556                                eps3dotdphij=eps3[0]*dphi[0][j]+eps3[1]*dphi[1][j]+eps3[2]*dphi[2][j];
     7557
     7558                                Ke->values[numdof*(4*i+0)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps1dotdphii;
     7559                                Ke->values[numdof*(4*i+0)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps1dotdphii;
     7560                                Ke->values[numdof*(4*i+0)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps1dotdphii;
     7561
     7562                                Ke->values[numdof*(4*i+1)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps2dotdphii;
     7563                                Ke->values[numdof*(4*i+1)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps2dotdphii;
     7564                                Ke->values[numdof*(4*i+1)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps2dotdphii;
     7565
     7566                                Ke->values[numdof*(4*i+2)+4*j+0]+=gauss->weight*Jdet*2*mu_prime*eps1dotdphij*eps3dotdphii;
     7567                                Ke->values[numdof*(4*i+2)+4*j+1]+=gauss->weight*Jdet*2*mu_prime*eps2dotdphij*eps3dotdphii;
     7568                                Ke->values[numdof*(4*i+2)+4*j+2]+=gauss->weight*Jdet*2*mu_prime*eps3dotdphij*eps3dotdphii;
     7569                        }
     7570                }
     7571        }
     7572
     7573        /*Transform Coordinate System*/
     7574        TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
     7575
     7576        /*Clean up and return*/
     7577        delete gauss;
     7578        return Ke;
     7579}
     7580/*}}}*/
    73217581/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHoriz{{{1*/
    73227582void  Penta::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
Note: See TracChangeset for help on using the changeset viewer.