Changeset 16982


Ignore:
Timestamp:
11/29/13 20:27:21 (11 years ago)
Author:
Mathieu Morlighem
Message:

CHG: use analysis to create P Vector instead of element

Location:
issm/trunk-jpl/src/c
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r16956 r16982  
    105105                virtual void       SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;
    106106                virtual ElementMatrix* CreateKMatrix(void)=0;
    107                 virtual ElementVector* CreatePVector(void)=0;
    108107                virtual void   CreateDVector(Vector<IssmDouble>* df)=0;
    109108                virtual void   CreateJacobianMatrix(Matrix<IssmDouble>* Jff)=0;

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

    r16974 r16982  
    540540                delete De;
    541541        }
    542 }
    543 /*}}}*/
    544 /*FUNCTION Penta::CreatePVector(void) {{{*/
    545 ElementVector* Penta::CreatePVector(void){
    546 
    547         /*retrieve parameters: */
    548         int analysis_type;
    549         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    550 
    551         /*if debugging mode, check that all pointers exist {{{*/
    552         _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs);
    553         /*}}}*/
    554 
    555         /*Skip if water element*/
    556         if(NoIceInElement()) return NULL;
    557 
    558         /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
    559         switch(analysis_type){
    560                 #ifdef _HAVE_STRESSBALANCE_
    561                 case StressbalanceAnalysisEnum:
    562                         return CreatePVectorStressbalanceHoriz();
    563                         break;
    564                 case StressbalanceSIAAnalysisEnum:
    565                         return CreatePVectorStressbalanceSIA();
    566                         break;
    567                 case StressbalanceVerticalAnalysisEnum:
    568                         return CreatePVectorStressbalanceVert();
    569                         break;
    570                 #endif
    571                 #ifdef _HAVE_CONTROL_
    572                 case AdjointHorizAnalysisEnum:
    573                         return CreatePVectorAdjointHoriz();
    574                         break;
    575                 #endif
    576                 #ifdef _HAVE_THERMAL_
    577                 case ThermalAnalysisEnum:
    578                         return CreatePVectorThermal();
    579                         break;
    580                 case EnthalpyAnalysisEnum:
    581                         return CreatePVectorEnthalpy();
    582                         break;
    583                 case MeltingAnalysisEnum:
    584                         return CreatePVectorMelting();
    585                         break;
    586                 #endif
    587                 case L2ProjectionBaseAnalysisEnum:
    588                         return CreatePVectorL2ProjectionBase();
    589                         break;
    590                 case MasstransportAnalysisEnum:
    591                         return CreatePVectorMasstransport();
    592                         break;
    593                 case FreeSurfaceTopAnalysisEnum:
    594                         return CreatePVectorFreeSurfaceTop();
    595                         break;
    596                 case FreeSurfaceBaseAnalysisEnum:
    597                         return CreatePVectorFreeSurfaceBase();
    598                         break;
    599                 #ifdef _HAVE_BALANCED_
    600                 case BalancethicknessAnalysisEnum:
    601                         return CreatePVectorBalancethickness();
    602                         break;
    603                 #endif
    604                 #ifdef _HAVE_HYDROLOGY_
    605           case HydrologyDCInefficientAnalysisEnum:
    606                         return CreatePVectorHydrologyDCInefficient();
    607                         break;
    608           case HydrologyDCEfficientAnalysisEnum:
    609                         return CreatePVectorHydrologyDCEfficient();
    610                         break;
    611           case L2ProjectionEPLAnalysisEnum:
    612                         return CreatePVectorL2ProjectionEPL();
    613                         break;
    614                 #endif
    615                 default:
    616                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    617         }
    618 
    619         /*make compiler happy*/
    620         return NULL;
    621 }
    622 /*}}}*/
    623 /*FUNCTION Penta::CreatePVectorMasstransport {{{*/
    624 ElementVector* Penta::CreatePVectorMasstransport(void){
    625 
    626         if (!IsOnBed()) return NULL;
    627 
    628         /*Depth Averaging Vx and Vy*/
    629         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    630         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    631 
    632         /*Call Tria function*/
    633         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    634         ElementVector* pe=tria->CreatePVectorMasstransport();
    635         delete tria->material; delete tria;
    636 
    637         /*Delete Vx and Vy averaged*/
    638         this->inputs->DeleteInput(VxAverageEnum);
    639         this->inputs->DeleteInput(VyAverageEnum);
    640 
    641         /*Clean up and return*/
    642         return pe;
    643 }
    644 /*}}}*/
    645 /*FUNCTION Penta::CreatePVectorFreeSurfaceTop {{{*/
    646 ElementVector* Penta::CreatePVectorFreeSurfaceTop(void){
    647 
    648         if(!IsOnSurface()) return NULL;
    649 
    650         /*Call Tria function*/
    651         Tria* tria=(Tria*)SpawnTria(1); //lower face is 0, upper face is 1.
    652         ElementVector* pe=tria->CreatePVectorFreeSurfaceTop();
    653         delete tria->material; delete tria;
    654 
    655         /*Clean up and return*/
    656         return pe;
    657 }
    658 /*}}}*/
    659 /*FUNCTION Penta::CreatePVectorFreeSurfaceBase {{{*/
    660 ElementVector* Penta::CreatePVectorFreeSurfaceBase(void){
    661 
    662         if(!IsOnBed()) return NULL;
    663 
    664         /*Call Tria function*/
    665         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    666         ElementVector* pe=tria->CreatePVectorFreeSurfaceBase();
    667         delete tria->material; delete tria;
    668 
    669         /*Clean up and return*/
    670         return pe;
    671 }
    672 /*}}}*/
    673 /*FUNCTION Penta::CreatePVectorL2ProjectionBase {{{*/
    674 ElementVector* Penta::CreatePVectorL2ProjectionBase(void){
    675 
    676         if (!IsOnBed()) return NULL;
    677 
    678         /*Call Tria function*/
    679         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    680         ElementVector* pe=tria->CreatePVectorL2Projection();
    681         delete tria->material; delete tria;
    682 
    683         /*clean up and return*/
    684         return pe;
    685542}
    686543/*}}}*/
    … …  
    40313888}
    40323889/*}}}*/
    4033 /*FUNCTION Penta::CreatePVectorEnthalpy {{{*/
    4034 ElementVector* Penta::CreatePVectorEnthalpy(void){
    4035 
    4036         /*compute all load vectors for this element*/
    4037         ElementVector* pe1=CreatePVectorEnthalpyVolume();
    4038         ElementVector* pe2=CreatePVectorEnthalpySheet();
    4039         ElementVector* pe3=CreatePVectorEnthalpyShelf();
    4040         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    4041 
    4042         /*clean-up and return*/
    4043         delete pe1;
    4044         delete pe2;
    4045         delete pe3;
    4046         return pe;
    4047 }
    4048 /*}}}*/
    4049 /*FUNCTION Penta::CreatePVectorEnthalpyVolume {{{*/
    4050 ElementVector* Penta::CreatePVectorEnthalpyVolume(void){
    4051 
    4052         /*Constants*/
    4053         const int  numdof=NUMVERTICES*NDOF1;
    4054 
    4055         /*Intermediaries*/
    4056         int        i;
    4057         int        stabilization;
    4058         IssmDouble Jdet,phi,dt;
    4059         IssmDouble rho_ice,heatcapacity;
    4060         IssmDouble thermalconductivity,kappa;
    4061         IssmDouble viscosity,pressure,enthalpy;
    4062         IssmDouble enthalpypicard[NUMVERTICES], pressure_p[NUMVERTICES];
    4063         IssmDouble tau_parameter,diameter;
    4064         IssmDouble u,v,w;
    4065         IssmDouble scalar_def,scalar_transient;
    4066         IssmDouble xyz_list[NUMVERTICES][3];
    4067         IssmDouble L[numdof];
    4068         IssmDouble dbasis[3][6];
    4069         IssmDouble epsilon[6];
    4070         GaussPenta *gauss=NULL;
    4071 
    4072         /*Initialize Element vector*/
    4073         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4074 
    4075         /*Retrieve all inputs and parameters*/
    4076         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4077         rho_ice=matpar->GetRhoIce();
    4078         heatcapacity=matpar->GetHeatCapacity();
    4079         thermalconductivity=matpar->GetThermalConductivity();
    4080         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4081         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    4082         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    4083         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    4084         Input* vz_input=inputs->GetInput(VzEnum);                                  _assert_(vz_input);
    4085         Input* enthalpy_input=NULL;
    4086         if(reCast<bool,IssmDouble>(dt)){
    4087                 enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input);
    4088         }
    4089         if (stabilization==2){
    4090                 diameter=MinEdgeLength(&xyz_list[0][0]);
    4091         }
    4092 
    4093         /* Start  looping on the number of gaussian points: */
    4094         gauss=new GaussPenta(3,3);
    4095         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4096 
    4097                 gauss->GaussPoint(ig);
    4098 
    4099                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    4100                 GetNodalFunctionsP1(&L[0], gauss);
    4101 
    4102                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    4103                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    4104                 GetPhi(&phi, &epsilon[0], viscosity);
    4105 
    4106                 scalar_def=phi/rho_ice*Jdet*gauss->weight;
    4107                 if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
    4108                
    4109                 /*TODO: add -beta*laplace T_m(p)*/
    4110                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
    4111 
    4112                 /* Build transient now */
    4113                 if(reCast<bool,IssmDouble>(dt)){
    4114                         enthalpy_input->GetInputValue(&enthalpy, gauss);
    4115                         scalar_transient=enthalpy*Jdet*gauss->weight;
    4116                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
    4117                 }
    4118 
    4119                 if(stabilization==2){
    4120                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    4121 
    4122                         vx_input->GetInputValue(&u, gauss);
    4123                         vy_input->GetInputValue(&v, gauss);
    4124                         vz_input->GetInputValue(&w, gauss);
    4125                         GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum);
    4126                         GetInputListOnVertices(&pressure_p[0],PressureEnum);
    4127                         kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpypicard,pressure_p);
    4128                         kappa/=rho_ice;
    4129                         tau_parameter=StabilizationParameter(u,v,w,diameter,kappa);
    4130 
    4131                         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]);
    4132                         if(reCast<bool,IssmDouble>(dt)){
    4133                                 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]);
    4134                         }
    4135                 }
    4136         }
    4137 
    4138         /*Clean up and return*/
    4139         delete gauss;
    4140         return pe;
    4141 }
    4142 /*}}}*/
    4143 /*FUNCTION Penta::CreatePVectorEnthalpyShelf {{{*/
    4144 ElementVector* Penta::CreatePVectorEnthalpyShelf(void){
    4145 
    4146         /*Constants*/
    4147         const int  numdof=NUMVERTICES*NDOF1;
    4148 
    4149         /*Intermediaries */
    4150         int        i,j;
    4151         IssmDouble Jdet2d;
    4152         IssmDouble heatcapacity,h_pmp;
    4153         IssmDouble mixed_layer_capacity,thermal_exchange_velocity;
    4154         IssmDouble rho_ice,rho_water,pressure,dt,scalar_ocean;
    4155         IssmDouble xyz_list[NUMVERTICES][3];
    4156         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    4157         IssmDouble basis[NUMVERTICES];
    4158         GaussPenta* gauss=NULL;
    4159 
    4160         /* Ice/ocean heat exchange flux on ice shelf base */
    4161         if (!IsOnBed() || !IsFloating()) return NULL;
    4162 
    4163         /*Initialize Element vector*/
    4164         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4165 
    4166         /*Retrieve all inputs and parameters*/
    4167         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4168         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4169         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    4170         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    4171         rho_water=matpar->GetRhoWater();
    4172         rho_ice=matpar->GetRhoIce();
    4173         heatcapacity=matpar->GetHeatCapacity();
    4174         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4175         Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    4176 
    4177         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    4178         gauss=new GaussPenta(0,1,2,2);
    4179         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4180 
    4181                 gauss->GaussPoint(ig);
    4182 
    4183                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    4184                 GetNodalFunctionsP1(&basis[0], gauss);
    4185 
    4186                 pressure_input->GetInputValue(&pressure,gauss);
    4187                 h_pmp=matpar->PureIceEnthalpy(pressure);
    4188 
    4189                 scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(h_pmp)/(rho_ice*heatcapacity);
    4190                 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
    4191 
    4192                 for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
    4193         }
    4194 
    4195         /*Clean up and return*/
    4196         delete gauss;
    4197         return pe;
    4198 }
    4199 /*}}}*/
    4200 /*FUNCTION Penta::CreatePVectorEnthalpySheet {{{*/
    4201 ElementVector* Penta::CreatePVectorEnthalpySheet(void){
    4202 
    4203         /*Constants*/
    4204         const int  numdof=NUMVERTICES*NDOF1;
    4205 
    4206         /*Intermediaries */
    4207         int        i,j;
    4208         int        analysis_type;
    4209         IssmDouble xyz_list[NUMVERTICES][3];
    4210         IssmDouble xyz_list_tria[NUMVERTICES2D][3]={0.0};
    4211         IssmDouble Jdet2d,dt;
    4212         IssmDouble rho_ice,heatcapacity;
    4213         IssmDouble geothermalflux_value, heatflux_value;
    4214         IssmDouble basalfriction,alpha2,vx,vy,vz;
    4215         IssmDouble scalar,enthalpy,enthalpyup;
    4216         IssmDouble pressure,pressureup;
    4217         IssmDouble watercolumn;
    4218         IssmDouble basis[NUMVERTICES];
    4219         Friction*  friction=NULL;
    4220         GaussPenta* gauss=NULL;
    4221         GaussPenta* gaussup=NULL;
    4222 
    4223         /* Geothermal flux on ice sheet base and basal friction */
    4224         if (!IsOnBed() || IsFloating()) return NULL;
    4225 
    4226         /*Initialize Element vector*/
    4227         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4228 
    4229         /*Retrieve all inputs and parameters*/
    4230         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4231         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4232         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4233         rho_ice=matpar->GetRhoIce();
    4234         heatcapacity=matpar->GetHeatCapacity();
    4235         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4236         Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
    4237         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
    4238         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
    4239         Input* enthalpy_input=inputs->GetInput(EnthalpyPicardEnum);       _assert_(enthalpy_input);
    4240         Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
    4241         Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    4242         Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);       _assert_(watercolumn_input);
    4243 
    4244         /*Build friction element, needed later: */
    4245         friction=new Friction(this,3);
    4246 
    4247         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    4248         gauss=new GaussPenta(0,1,2,2);
    4249         gaussup=new GaussPenta(3,4,5,2);
    4250 
    4251         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4252 
    4253                 gauss->GaussPoint(ig);
    4254                 gaussup->GaussPoint(ig);
    4255 
    4256                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    4257                 GetNodalFunctionsP1(&basis[0], gauss);
    4258 
    4259                 enthalpy_input->GetInputValue(&enthalpy,gauss);
    4260                 pressure_input->GetInputValue(&pressure,gauss);
    4261                 watercolumn_input->GetInputValue(&watercolumn,gauss);
    4262 
    4263                 if((watercolumn<=0.) && (enthalpy<matpar->PureIceEnthalpy(pressure))){
    4264                         /* the above check is equivalent to
    4265                          NOT ((watercolumn>0.) AND (enthalpy<PIE)) AND (enthalpy<PIE)*/
    4266                         geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4267 
    4268                         friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    4269                         vx_input->GetInputValue(&vx,gauss);
    4270                         vy_input->GetInputValue(&vy,gauss);
    4271                         vz_input->GetInputValue(&vz,gauss);
    4272                         basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)+pow(vz,2.0));
    4273 
    4274                         heatflux_value=(basalfriction+geothermalflux_value)/(rho_ice);
    4275                         scalar=gauss->weight*Jdet2d*heatflux_value;
    4276                         if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
    4277                         for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    4278                 }
    4279                 else if(enthalpy>=matpar->PureIceEnthalpy(pressure)){
    4280                         /* check positive thickness of temperate basal ice layer */
    4281                         enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    4282                         pressure_input->GetInputValue(&pressureup,gaussup);
    4283                         if(enthalpyup>=matpar->PureIceEnthalpy(pressureup)){
    4284                                 // TODO: temperate ice has positive thickness: grad enthalpy*n=0.
    4285                         }
    4286                         else{
    4287                                 // only base temperate, set Dirichlet BCs in Penta::UpdateBasalConstraintsEnthalpy()
    4288                         }
    4289                 }
    4290                 else {
    4291                         // base cold, but watercolumn positive. Set base to h_pmp.
    4292                 }
    4293         }
    4294 
    4295         /*Clean up and return*/
    4296         delete gauss;
    4297         delete gaussup;
    4298         delete friction;
    4299         return pe;
    4300 }
    4301 /*}}}*/
    4302 /*FUNCTION Penta::CreatePVectorMelting {{{*/
    4303 ElementVector* Penta::CreatePVectorMelting(void){
    4304         return NULL;
    4305 }
    4306 /*}}}*/
    4307 /*FUNCTION Penta::CreatePVectorThermal {{{*/
    4308 ElementVector* Penta::CreatePVectorThermal(void){
    4309 
    4310         /*compute all load vectors for this element*/
    4311         ElementVector* pe1=CreatePVectorThermalVolume();
    4312         ElementVector* pe2=CreatePVectorThermalSheet();
    4313         ElementVector* pe3=CreatePVectorThermalShelf();
    4314         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    4315 
    4316         /*clean-up and return*/
    4317         delete pe1;
    4318         delete pe2;
    4319         delete pe3;
    4320         return pe;
    4321 }
    4322 /*}}}*/
    4323 /*FUNCTION Penta::CreatePVectorThermalVolume {{{*/
    4324 ElementVector* Penta::CreatePVectorThermalVolume(void){
    4325 
    4326         /*Constants*/
    4327         const int  numdof=NUMVERTICES*NDOF1;
    4328 
    4329         /*Intermediaries*/
    4330         int        i;
    4331         int        stabilization;
    4332         IssmDouble Jdet,phi,dt;
    4333         IssmDouble rho_ice,heatcapacity;
    4334         IssmDouble thermalconductivity,kappa;
    4335         IssmDouble viscosity,temperature;
    4336         IssmDouble tau_parameter,diameter;
    4337         IssmDouble u,v,w;
    4338         IssmDouble scalar_def,scalar_transient;
    4339         IssmDouble xyz_list[NUMVERTICES][3];
    4340         IssmDouble L[numdof];
    4341         IssmDouble dbasis[3][6];
    4342         IssmDouble epsilon[6];
    4343         GaussPenta *gauss=NULL;
    4344 
    4345         /*Initialize Element vector*/
    4346         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4347 
    4348         /*Retrieve all inputs and parameters*/
    4349         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4350         rho_ice=matpar->GetRhoIce();
    4351         heatcapacity=matpar->GetHeatCapacity();
    4352         thermalconductivity=matpar->GetThermalConductivity();
    4353         kappa=thermalconductivity/(rho_ice*heatcapacity);
    4354         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4355         this->parameters->FindParam(&stabilization,ThermalStabilizationEnum);
    4356         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    4357         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    4358         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    4359         Input* temperature_input=NULL;
    4360         if(reCast<bool,IssmDouble>(dt)) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs);
    4361         if(stabilization==2) diameter=MinEdgeLength(&xyz_list[0][0]);
    4362 
    4363         /* Start  looping on the number of gaussian points: */
    4364         gauss=new GaussPenta(3,3);
    4365         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4366 
    4367                 gauss->GaussPoint(ig);
    4368 
    4369                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    4370                 GetNodalFunctionsP1(&L[0], gauss);
    4371 
    4372                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    4373                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    4374                 GetPhi(&phi, &epsilon[0], viscosity);
    4375 
    4376                 scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight;
    4377                 if(reCast<bool,IssmDouble>(dt)) scalar_def=scalar_def*dt;
    4378 
    4379                 for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_def*L[i];
    4380 
    4381                 /* Build transient now */
    4382                 if(reCast<bool,IssmDouble>(dt)){
    4383                         temperature_input->GetInputValue(&temperature, gauss);
    4384                         scalar_transient=temperature*Jdet*gauss->weight;
    4385                         for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=scalar_transient*L[i];
    4386                 }
    4387 
    4388                 if(stabilization==2){
    4389                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    4390 
    4391                         vx_input->GetInputValue(&u, gauss);
    4392                         vy_input->GetInputValue(&v, gauss);
    4393                         vz_input->GetInputValue(&w, gauss);
    4394 
    4395                         tau_parameter=StabilizationParameter(u,v,w,diameter,kappa);
    4396 
    4397                         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]);
    4398                         if(reCast<bool,IssmDouble>(dt)){
    4399                                 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]);
    4400                         }
    4401                 }
    4402         }
    4403 
    4404         /*Clean up and return*/
    4405         delete gauss;
    4406         return pe;
    4407 }
    4408 /*}}}*/
    4409 /*FUNCTION Penta::CreatePVectorThermalShelf {{{*/
    4410 ElementVector* Penta::CreatePVectorThermalShelf(void){
    4411 
    4412         /*Constants*/
    4413         const int  numdof=NUMVERTICES*NDOF1;
    4414 
    4415         /*Intermediaries */
    4416         int        i,j;
    4417         IssmDouble Jdet2d;
    4418         IssmDouble mixed_layer_capacity,thermal_exchange_velocity;
    4419         IssmDouble rho_ice,rho_water,pressure,dt,scalar_ocean;
    4420         IssmDouble heatcapacity,t_pmp;
    4421         IssmDouble xyz_list[NUMVERTICES][3];
    4422         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    4423         IssmDouble basis[NUMVERTICES];
    4424         GaussPenta* gauss=NULL;
    4425 
    4426         /* Ice/ocean heat exchange flux on ice shelf base */
    4427         if (!IsOnBed() || !IsFloating()) return NULL;
    4428 
    4429         /*Initialize Element vector*/
    4430         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4431 
    4432         /*Retrieve all inputs and parameters*/
    4433         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4434         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4435         mixed_layer_capacity=matpar->GetMixedLayerCapacity();
    4436         thermal_exchange_velocity=matpar->GetThermalExchangeVelocity();
    4437         rho_water=matpar->GetRhoWater();
    4438         rho_ice=matpar->GetRhoIce();
    4439         heatcapacity=matpar->GetHeatCapacity();
    4440         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4441         Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
    4442 
    4443         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    4444         gauss=new GaussPenta(0,1,2,2);
    4445         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4446 
    4447                 gauss->GaussPoint(ig);
    4448 
    4449                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    4450                 GetNodalFunctionsP1(&basis[0], gauss);
    4451 
    4452                 pressure_input->GetInputValue(&pressure,gauss);
    4453                 t_pmp=matpar->TMeltingPoint(pressure);
    4454 
    4455                 scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(t_pmp)/(heatcapacity*rho_ice);
    4456                 if(reCast<bool,IssmDouble>(dt)) scalar_ocean=dt*scalar_ocean;
    4457 
    4458                 for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
    4459         }
    4460 
    4461         /*Clean up and return*/
    4462         delete gauss;
    4463         return pe;
    4464 }
    4465 /*}}}*/
    4466 /*FUNCTION Penta::CreatePVectorThermalSheet {{{*/
    4467 ElementVector* Penta::CreatePVectorThermalSheet(void){
    4468 
    4469         /*Constants*/
    4470         const int  numdof=NUMVERTICES*NDOF1;
    4471 
    4472         /*Intermediaries */
    4473         int         i,j;
    4474         int         analysis_type;
    4475         IssmDouble  xyz_list[NUMVERTICES][3];
    4476         IssmDouble  xyz_list_tria[NUMVERTICES2D][3]={0.0};
    4477         IssmDouble  Jdet2d,dt;
    4478         IssmDouble  rho_ice,heatcapacity,geothermalflux_value;
    4479         IssmDouble  basalfriction,alpha2,vx,vy;
    4480         IssmDouble  basis[NUMVERTICES];
    4481         IssmDouble  scalar;
    4482         Friction*   friction=NULL;
    4483         GaussPenta* gauss=NULL;
    4484 
    4485         /* Geothermal flux on ice sheet base and basal friction */
    4486         if (!IsOnBed() || IsFloating()) return NULL;
    4487 
    4488         /*Initialize Element vector*/
    4489         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    4490 
    4491         /*Retrieve all inputs and parameters*/
    4492         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4493         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    4494         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    4495         rho_ice=matpar->GetRhoIce();
    4496         heatcapacity=matpar->GetHeatCapacity();
    4497         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    4498         Input* vx_input=inputs->GetInput(VxEnum);                         _assert_(vx_input);
    4499         Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
    4500         Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
    4501         Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    4502 
    4503         /*Build frictoin element, needed later: */
    4504         friction=new Friction(this,3);
    4505 
    4506         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    4507         gauss=new GaussPenta(0,1,2,2);
    4508         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4509 
    4510                 gauss->GaussPoint(ig);
    4511 
    4512                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    4513                 GetNodalFunctionsP1(&basis[0], gauss);
    4514 
    4515                         geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4516                         friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
    4517                         vx_input->GetInputValue(&vx,gauss);
    4518                         vy_input->GetInputValue(&vy,gauss);
    4519                         basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    4520 
    4521                         scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    4522                         if(reCast<bool,IssmDouble>(dt)) scalar=dt*scalar;
    4523 
    4524                         for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
    4525         }
    4526 
    4527         /*Clean up and return*/
    4528         delete gauss;
    4529         delete friction;
    4530         return pe;
    4531 }
    4532 /*}}}*/
    45333890/*FUNCTION Penta::UpdateBasalConstraintsEnthalpy{{{*/
    45343891void  Penta::UpdateBasalConstraintsEnthalpy(void){
    … …  
    51184475        delete gauss;
    51194476        return Ke;
    5120 }
    5121 /*}}}*/
    5122 /*FUNCTION Penta::CreatePVectorAdjointHoriz{{{*/
    5123 ElementVector* Penta::CreatePVectorAdjointHoriz(void){
    5124 
    5125         int approximation;
    5126         inputs->GetInputValue(&approximation,ApproximationEnum);
    5127 
    5128         switch(approximation){
    5129                 case SSAApproximationEnum:
    5130                         return CreatePVectorAdjointSSA();
    5131                 case HOApproximationEnum:
    5132                         return CreatePVectorAdjointHO();
    5133                 case NoneApproximationEnum:
    5134                         return NULL;
    5135                 case FSApproximationEnum:
    5136                         return CreatePVectorAdjointFS();
    5137                 default:
    5138                         _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    5139         }
    5140 }
    5141 /*}}}*/
    5142 /*FUNCTION Penta::CreatePVectorAdjointSSA{{{*/
    5143 ElementVector* Penta::CreatePVectorAdjointSSA(){
    5144 
    5145         if (!IsOnBed()) return NULL;
    5146 
    5147         /*Call Tria function*/
    5148         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    5149         ElementVector* pe=tria->CreatePVectorAdjointHoriz();
    5150         delete tria->material; delete tria;
    5151 
    5152         /*clean up and return*/
    5153         return pe;
    5154 }
    5155 /*}}}*/
    5156 /*FUNCTION Penta::CreatePVectorAdjointHO{{{*/
    5157 ElementVector* Penta::CreatePVectorAdjointHO(void){
    5158 
    5159         /*Nothing to be done if not on surface*/
    5160         if (!IsOnSurface()) return NULL;
    5161 
    5162         /*Intermediaries */
    5163         int        i,j,resp;
    5164         int       *responses=NULL;
    5165         int        num_responses;
    5166         IssmDouble Jdet2d;
    5167         IssmDouble obs_velocity_mag,velocity_mag;
    5168         IssmDouble dux,duy;
    5169         IssmDouble epsvel=2.220446049250313e-16;
    5170         IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
    5171         IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.;
    5172         IssmDouble vx,vy,vxobs,vyobs,weight;
    5173         IssmDouble xyz_list[NUMVERTICES][3];
    5174         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    5175 
    5176         /*Fetch number of nodes and dof for this finite element*/
    5177         int numnodes = this->NumberofNodes();
    5178 
    5179         /*Initialize Element vector*/
    5180         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    5181         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    5182 
    5183         /*Retrieve all inputs and parameters*/
    5184         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5185         this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    5186         this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    5187         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    5188         Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    5189         Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
    5190         Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
    5191         Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
    5192 
    5193         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    5194 
    5195         /*Get Surface if required by one response*/
    5196         for(resp=0;resp<num_responses;resp++){
    5197                 if(responses[resp]==SurfaceAverageVelMisfitEnum){
    5198                         inputs->GetInputValue(&S,SurfaceAreaEnum); break;
    5199                 }
    5200         }
    5201 
    5202         /* Start  looping on the number of gaussian points: */
    5203         GaussPenta* gauss=new GaussPenta(3,4,5,4);
    5204         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5205 
    5206                 gauss->GaussPoint(ig);
    5207 
    5208                 /* Get Jacobian determinant: */
    5209                 GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss);
    5210 
    5211                 /*Get all parameters at gaussian point*/
    5212                 vx_input->GetInputValue(&vx,gauss);
    5213                 vy_input->GetInputValue(&vy,gauss);
    5214                 vxobs_input->GetInputValue(&vxobs,gauss);
    5215                 vyobs_input->GetInputValue(&vyobs,gauss);
    5216                 GetNodalFunctions(basis, gauss);
    5217 
    5218                 /*Loop over all requested responses*/
    5219                 for(resp=0;resp<num_responses;resp++){
    5220 
    5221                         weights_input->GetInputValue(&weight,gauss,responses[resp]);
    5222 
    5223                         switch(responses[resp]){
    5224                                 case SurfaceAbsVelMisfitEnum:
    5225                                         /*
    5226                                          *      1  [           2              2 ]
    5227                                          * J = --- | (u - u   )  +  (v - v   )  |
    5228                                          *      2  [       obs            obs   ]
    5229                                          *
    5230                                          *        dJ
    5231                                          * DU = - -- = (u   - u )
    5232                                          *        du     obs
    5233                                          */
    5234                                         for(i=0;i<numnodes;i++){
    5235                                                 dux=vxobs-vx;
    5236                                                 duy=vyobs-vy;
    5237                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i];
    5238                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i];
    5239                                         }
    5240                                         break;
    5241                                 case SurfaceRelVelMisfitEnum:
    5242                                         /*
    5243                                          *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    5244                                          * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    5245                                          *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    5246                                          *              obs                        obs                     
    5247                                          *
    5248                                          *        dJ     \bar{v}^2
    5249                                          * DU = - -- = ------------- (u   - u )
    5250                                          *        du   (u   + eps)^2    obs
    5251                                          *               obs
    5252                                          */
    5253                                         for(i=0;i<numnodes;i++){
    5254                                                 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
    5255                                                 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
    5256                                                 dux=scalex*(vxobs-vx);
    5257                                                 duy=scaley*(vyobs-vy);
    5258                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i];
    5259                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i];
    5260                                         }
    5261                                         break;
    5262                                 case SurfaceLogVelMisfitEnum:
    5263                                         /*
    5264                                          *                 [        vel + eps     ] 2
    5265                                          * J = 4 \bar{v}^2 | log ( -----------  ) | 
    5266                                          *                 [       vel   + eps    ]
    5267                                          *                            obs
    5268                                          *
    5269                                          *        dJ                 2 * log(...)
    5270                                          * DU = - -- = - 4 \bar{v}^2 -------------  u
    5271                                          *        du                 vel^2 + eps
    5272                                          *           
    5273                                          */
    5274                                         for(i=0;i<numnodes;i++){
    5275                                                 velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
    5276                                                 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
    5277                                                 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    5278                                                 dux=scale*vx;
    5279                                                 duy=scale*vy;
    5280                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i];
    5281                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i];
    5282                                         }
    5283                                         break;
    5284                                 case SurfaceAverageVelMisfitEnum:
    5285                                         /*
    5286                                          *      1                    2              2
    5287                                          * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    5288                                          *      S                obs            obs
    5289                                          *
    5290                                          *        dJ      1       1
    5291                                          * DU = - -- = - --- ----------- * 2 (u - u   )
    5292                                          *        du      S  2 sqrt(...)           obs
    5293                                          */
    5294                                         for(i=0;i<numnodes;i++){
    5295                                                 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
    5296                                                 dux=scale*(vxobs-vx);
    5297                                                 duy=scale*(vyobs-vy);
    5298                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i];
    5299                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i];
    5300                                         }
    5301                                         break;
    5302                                 case SurfaceLogVxVyMisfitEnum:
    5303                                         /*
    5304                                          *      1            [        |u| + eps     2          |v| + eps     2  ]
    5305                                          * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    5306                                          *      2            [       |u    |+ eps              |v    |+ eps     ]
    5307                                          *                              obs                       obs
    5308                                          *        dJ                              1      u                             1
    5309                                          * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    5310                                          *        du                         |u| + eps  |u|                           u + eps
    5311                                          */
    5312                                         for(i=0;i<numnodes;i++){
    5313                                                 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    5314                                                 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    5315                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet2d*gauss->weight*basis[i];
    5316                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet2d*gauss->weight*basis[i];
    5317                                         }
    5318                                         break;
    5319                                 case DragCoefficientAbsGradientEnum:
    5320                                         /*Nothing in P vector*/
    5321                                         break;
    5322                                 case ThicknessAbsGradientEnum:
    5323                                         /*Nothing in P vector*/
    5324                                         break;
    5325                                 case ThicknessAlongGradientEnum:
    5326                                         /*Nothing in P vector*/
    5327                                         break;
    5328                                 case ThicknessAcrossGradientEnum:
    5329                                         /*Nothing in P vector*/
    5330                                         break;
    5331                                 case RheologyBbarAbsGradientEnum:
    5332                                         /*Nothing in P vector*/
    5333                                         break;
    5334                                 default:
    5335                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    5336                         }
    5337                 }
    5338         }
    5339 
    5340         /*Clean up and return*/
    5341         xDelete<IssmDouble>(basis);
    5342         xDelete<int>(responses);
    5343         delete gauss;
    5344         return pe;
    5345 }
    5346 /*}}}*/
    5347 /*FUNCTION Penta::CreatePVectorAdjointFS{{{*/
    5348 ElementVector* Penta::CreatePVectorAdjointFS(void){
    5349 
    5350         if(!IsOnSurface()) return NULL;
    5351 
    5352         /*Intermediaries */
    5353         int         i,j,resp;
    5354         int        *responses= NULL;
    5355         int         num_responses;
    5356         IssmDouble  Jdet2d;
    5357         IssmDouble  obs_velocity_mag,velocity_mag;
    5358         IssmDouble  dux,duy;
    5359         IssmDouble  epsvel  = 2.220446049250313e-16;
    5360         IssmDouble  meanvel = 3.170979198376458e-05;  /*1000 m/yr */
    5361         IssmDouble  scalex  = 0,scaley=0,scale=0,S=0;
    5362         IssmDouble  vx,vy,vxobs,vyobs,weight;
    5363         IssmDouble  xyz_list[NUMVERTICES][3];
    5364         IssmDouble      xyz_list_tria[NUMVERTICES2D][3];
    5365 
    5366         /*Fetch number of nodes and dof for this finite element*/
    5367         int vnumnodes = this->NumberofNodesVelocity();
    5368         int pnumnodes = this->NumberofNodesPressure();
    5369 
    5370         /*Prepare coordinate system list*/
    5371         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    5372         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    5373         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    5374 
    5375         /*Initialize Element matrix and vectors*/
    5376         ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
    5377         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    5378 
    5379         /*Retrieve all inputs and parameters*/
    5380         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5381         this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    5382         this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    5383         Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    5384         Input* vx_input      = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
    5385         Input* vy_input      = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
    5386         Input* vxobs_input   = inputs->GetInput(InversionVxObsEnum);                     _assert_(vxobs_input);
    5387         Input* vyobs_input   = inputs->GetInput(InversionVyObsEnum);                     _assert_(vyobs_input);
    5388 
    5389         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    5390 
    5391         /*Get Surface if required by one response*/
    5392         for(resp=0;resp<num_responses;resp++){
    5393                 if(responses[resp]==SurfaceAverageVelMisfitEnum){
    5394                         inputs->GetInputValue(&S,SurfaceAreaEnum); break;
    5395                 }
    5396         }
    5397 
    5398         /* Start  looping on the number of gaussian points: */
    5399         GaussPenta* gauss=new GaussPenta(3,4,5,4);
    5400         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5401 
    5402                 gauss->GaussPoint(ig);
    5403 
    5404                 /* Get Jacobian determinant: */
    5405                 GetTriaJacobianDeterminant(&Jdet2d,&xyz_list_tria[0][0],gauss);
    5406 
    5407                 /*Get all parameters at gaussian point*/
    5408                 vx_input->GetInputValue(&vx,gauss);
    5409                 vy_input->GetInputValue(&vy,gauss);
    5410                 vxobs_input->GetInputValue(&vxobs,gauss);
    5411                 vyobs_input->GetInputValue(&vyobs,gauss);
    5412                 GetNodalFunctionsVelocity(vbasis,gauss);
    5413 
    5414                 /*Loop over all requested responses*/
    5415                 for(resp=0;resp<num_responses;resp++){
    5416 
    5417                         weights_input->GetInputValue(&weight,gauss,responses[resp]);
    5418 
    5419                         switch(responses[resp]){
    5420 
    5421                                 case SurfaceAbsVelMisfitEnum:
    5422                                         /*
    5423                                          *      1  [           2              2 ]
    5424                                          * J = --- | (u - u   )  +  (v - v   )  |
    5425                                          *      2  [       obs            obs   ]
    5426                                          *
    5427                                          *        dJ
    5428                                          * DU = - -- = (u   - u )
    5429                                          *        du     obs
    5430                                          */
    5431                                         for(i=0;i<vnumnodes;i++){
    5432                                                 dux=vxobs-vx;
    5433                                                 duy=vyobs-vy;
    5434                                                 pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
    5435                                                 pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
    5436                                         }
    5437                                         break;
    5438                                 case SurfaceRelVelMisfitEnum:
    5439                                         /*
    5440                                          *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    5441                                          * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    5442                                          *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    5443                                          *              obs                        obs                     
    5444                                          *
    5445                                          *        dJ     \bar{v}^2
    5446                                          * DU = - -- = ------------- (u   - u )
    5447                                          *        du   (u   + eps)^2    obs
    5448                                          *               obs
    5449                                          */
    5450                                         for(i=0;i<vnumnodes;i++){
    5451                                                 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
    5452                                                 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
    5453                                                 dux=scalex*(vxobs-vx);
    5454                                                 duy=scaley*(vyobs-vy);
    5455                                                 pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
    5456                                                 pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
    5457                                         }
    5458                                         break;
    5459                                 case SurfaceLogVelMisfitEnum:
    5460                                         /*
    5461                                          *                 [        vel + eps     ] 2
    5462                                          * J = 4 \bar{v}^2 | log ( -----------  ) | 
    5463                                          *                 [       vel   + eps    ]
    5464                                          *                            obs
    5465                                          *
    5466                                          *        dJ                 2 * log(...)
    5467                                          * DU = - -- = - 4 \bar{v}^2 -------------  u
    5468                                          *        du                 vel^2 + eps
    5469                                          *           
    5470                                          */
    5471                                         for(i=0;i<vnumnodes;i++){
    5472                                                 velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
    5473                                                 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
    5474                                                 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    5475                                                 dux=scale*vx;
    5476                                                 duy=scale*vy;
    5477                                                 pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
    5478                                                 pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
    5479                                         }
    5480                                         break;
    5481                                 case SurfaceAverageVelMisfitEnum:
    5482                                         /*
    5483                                          *      1                    2              2
    5484                                          * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    5485                                          *      S                obs            obs
    5486                                          *
    5487                                          *        dJ      1       1
    5488                                          * DU = - -- = - --- ----------- * 2 (u - u   )
    5489                                          *        du      S  2 sqrt(...)           obs
    5490                                          */
    5491                                         for(i=0;i<vnumnodes;i++){
    5492                                                 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
    5493                                                 dux=scale*(vxobs-vx);
    5494                                                 duy=scale*(vyobs-vy);
    5495                                                 pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
    5496                                                 pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
    5497                                         }
    5498                                         break;
    5499                                 case SurfaceLogVxVyMisfitEnum:
    5500                                         /*
    5501                                          *      1            [        |u| + eps     2          |v| + eps     2  ]
    5502                                          * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    5503                                          *      2            [       |u    |+ eps              |v    |+ eps     ]
    5504                                          *                              obs                       obs
    5505                                          *        dJ                              1      u                             1
    5506                                          * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    5507                                          *        du                         |u| + eps  |u|                           u + eps
    5508                                          */
    5509                                         for(i=0;i<vnumnodes;i++){
    5510                                                 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    5511                                                 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    5512                                                 pe->values[i*NDOF3+0]+=dux*weight*Jdet2d*gauss->weight*vbasis[i];
    5513                                                 pe->values[i*NDOF3+1]+=duy*weight*Jdet2d*gauss->weight*vbasis[i];
    5514                                         }
    5515                                         break;
    5516                                 case DragCoefficientAbsGradientEnum:
    5517                                         /*Nothing in P vector*/
    5518                                         break;
    5519                                 case ThicknessAbsGradientEnum:
    5520                                         /*Nothing in P vector*/
    5521                                         break;
    5522                                 case ThicknessAcrossGradientEnum:
    5523                                         /*Nothing in P vector*/
    5524                                         break;
    5525                                 case ThicknessAlongGradientEnum:
    5526                                         /*Nothing in P vector*/
    5527                                         break;
    5528                                 case RheologyBbarAbsGradientEnum:
    5529                                         /*Nothing in P vector*/
    5530                                         break;
    5531                                 default:
    5532                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    5533                         }
    5534                 }
    5535         }
    5536 
    5537         /*Clean up and return*/
    5538         xDelete<int>(responses);
    5539         xDelete<int>(cs_list);
    5540         xDelete<IssmDouble>(vbasis);
    5541         delete gauss;
    5542         return pe;
    55434477}
    55444478/*}}}*/
    … …  
    78676801}
    78686802/*}}}*/
    7869 /*FUNCTION Penta::CreatePVectorCouplingSSAFS {{{*/
    7870 ElementVector* Penta::CreatePVectorCouplingSSAFS(void){
    7871 
    7872         /*compute all load vectors for this element*/
    7873         ElementVector* pe1=CreatePVectorCouplingSSAFSViscous();
    7874         ElementVector* pe2=CreatePVectorCouplingSSAFSFriction();
    7875         ElementVector* pe =new ElementVector(pe1,pe2);
    7876 
    7877         /*clean-up and return*/
    7878         delete pe1;
    7879         delete pe2;
    7880         return pe;
    7881 }
    7882 /*}}}*/
    7883 /*FUNCTION Penta::CreatePVectorCouplingSSAFSViscous {{{*/
    7884 ElementVector* Penta::CreatePVectorCouplingSSAFSViscous(void){
    7885 
    7886         /*Intermediaries */
    7887         int         i,approximation;
    7888         IssmDouble  viscosity,Jdet;
    7889         IssmDouble  FSreconditioning;
    7890         IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    7891         IssmDouble  dw[3];
    7892         IssmDouble  xyz_list[NUMVERTICES][3];
    7893         IssmDouble  basis[6]; //for the six nodes of the penta
    7894         IssmDouble  dbasis[3][6]; //for the six nodes of the penta
    7895         GaussPenta *gauss=NULL;
    7896 
    7897         /*Initialize Element vector and return if necessary*/
    7898         inputs->GetInputValue(&approximation,ApproximationEnum);
    7899         if(approximation!=SSAFSApproximationEnum) return NULL;
    7900         int vnumnodes = this->NumberofNodesVelocity();
    7901         int pnumnodes = this->NumberofNodesPressure();
    7902 
    7903         /*Prepare coordinate system list*/
    7904         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    7905         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    7906         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    7907         ElementVector* pe=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    7908 
    7909         /*Retrieve all inputs and parameters*/
    7910         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7911         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    7912         Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    7913         Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    7914         Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    7915         Input* vzSSA_input=inputs->GetInput(VzSSAEnum);   _assert_(vzSSA_input);
    7916 
    7917         /* Start  looping on the number of gaussian points: */
    7918         gauss=new GaussPenta(5,5);
    7919         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7920 
    7921                 gauss->GaussPoint(ig);
    7922 
    7923                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7924                 GetNodalFunctionsP1(&basis[0], gauss);
    7925                 GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    7926 
    7927                 vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    7928 
    7929                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    7930                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    7931 
    7932                 for(i=0;i<NUMVERTICES;i++){
    7933                         pe->values[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
    7934                         pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    7935                         pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    7936                         pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    7937                 }
    7938         }
    7939 
    7940         /*Transform coordinate system*/
    7941         ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    7942 
    7943         /*Clean up and return*/
    7944         xDelete<int>(cs_list);
    7945         delete gauss;
    7946         return pe;
    7947 }
    7948 /*}}}*/
    7949 /*FUNCTION Penta::CreatePVectorCouplingSSAFSFriction{{{*/
    7950 ElementVector* Penta::CreatePVectorCouplingSSAFSFriction(void){
    7951 
    7952         /*Intermediaries*/
    7953         int         i,j;
    7954         int         approximation,analysis_type;
    7955         IssmDouble  Jdet,Jdet2d;
    7956         IssmDouble  FSreconditioning;
    7957         IssmDouble      bed_normal[3];
    7958         IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    7959         IssmDouble  viscosity, w, alpha2_gauss;
    7960         IssmDouble  dw[3];
    7961         IssmDouble      xyz_list_tria[NUMVERTICES2D][3];
    7962         IssmDouble  xyz_list[NUMVERTICES][3];
    7963         IssmDouble  basis[6]; //for the six nodes of the penta
    7964         Friction*   friction=NULL;
    7965         GaussPenta  *gauss=NULL;
    7966 
    7967         /*Initialize Element vector and return if necessary*/
    7968         if(!IsOnBed() || IsFloating()) return NULL;
    7969         inputs->GetInputValue(&approximation,ApproximationEnum);
    7970         if(approximation!=SSAFSApproximationEnum) return NULL;
    7971         int vnumnodes = this->NumberofNodesVelocity();
    7972         int pnumnodes = this->NumberofNodesPressure();
    7973 
    7974         /*Prepare coordinate system list*/
    7975         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    7976         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    7977         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    7978         ElementVector* pe=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    7979 
    7980         /*Retrieve all inputs and parameters*/
    7981         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7982         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    7983         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    7984         Input* vx_input=inputs->GetInput(VxEnum);        _assert_(vx_input);
    7985         Input* vy_input=inputs->GetInput(VyEnum);        _assert_(vy_input);
    7986         Input* vz_input=inputs->GetInput(VzEnum);        _assert_(vz_input);
    7987         Input* vzSSA_input=inputs->GetInput(VzSSAEnum);  _assert_(vzSSA_input);
    7988 
    7989         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    7990 
    7991         /*build friction object, used later on: */
    7992         friction=new Friction(this,3);
    7993 
    7994         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    7995         gauss=new GaussPenta(0,1,2,2);
    7996         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7997 
    7998                 gauss->GaussPoint(ig);
    7999 
    8000                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    8001                 GetNodalFunctionsP1(basis, gauss);
    8002 
    8003                 vzSSA_input->GetInputValue(&w, gauss);
    8004                 vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    8005 
    8006                 NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    8007                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8008                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    8009                 friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    8010 
    8011                 for(i=0;i<NUMVERTICES2D;i++){
    8012                         pe->values[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
    8013                         pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
    8014                         pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
    8015                 }
    8016         }
    8017 
    8018         /*Transform coordinate system*/
    8019         ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    8020 
    8021         /*Clean up and return*/
    8022         xDelete<int>(cs_list);
    8023         delete gauss;
    8024         delete friction;
    8025         return pe;
    8026 }
    8027 /*}}}*/
    8028 /*FUNCTION Penta::CreatePVectorCouplingHOFS {{{*/
    8029 ElementVector* Penta::CreatePVectorCouplingHOFS(void){
    8030 
    8031         /*compute all load vectors for this element*/
    8032         ElementVector* pe1=CreatePVectorCouplingHOFSViscous();
    8033         ElementVector* pe2=CreatePVectorCouplingHOFSFriction();
    8034         ElementVector* pe =new ElementVector(pe1,pe2);
    8035 
    8036         /*clean-up and return*/
    8037         delete pe1;
    8038         delete pe2;
    8039         return pe;
    8040 }
    8041 /*}}}*/
    8042 /*FUNCTION Penta::CreatePVectorCouplingHOFSViscous {{{*/
    8043 ElementVector* Penta::CreatePVectorCouplingHOFSViscous(void){
    8044 
    8045         /*Intermediaries */
    8046         int         i;
    8047         int         approximation;
    8048         IssmDouble  viscosity,Jdet;
    8049         IssmDouble  FSreconditioning;
    8050         IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    8051         IssmDouble  dw[3];
    8052         IssmDouble  xyz_list[NUMVERTICES][3];
    8053         IssmDouble  basis[6]; //for the six nodes of the penta
    8054         IssmDouble  dbasis[3][6]; //for the six nodes of the penta
    8055         GaussPenta *gauss=NULL;
    8056 
    8057         /*Initialize Element vector and return if necessary*/
    8058         inputs->GetInputValue(&approximation,ApproximationEnum);
    8059         if(approximation!=HOFSApproximationEnum) return NULL;
    8060         int vnumnodes = this->NumberofNodesVelocity();
    8061         int pnumnodes = this->NumberofNodesPressure();
    8062 
    8063         /*Prepare coordinate system list*/
    8064         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    8065         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    8066         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    8067 
    8068         /*Initialize Element matrix and vectors*/
    8069         ElementVector* pe=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    8070 
    8071         /*Retrieve all inputs and parameters*/
    8072         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8073         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    8074         Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
    8075         Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
    8076         Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
    8077         Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
    8078 
    8079         /* Start  looping on the number of gaussian points: */
    8080         gauss=new GaussPenta(5,5);
    8081         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8082 
    8083                 gauss->GaussPoint(ig);
    8084 
    8085                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8086                 GetNodalFunctionsP1(&basis[0], gauss);
    8087                 GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
    8088 
    8089                 vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    8090 
    8091                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8092                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    8093 
    8094                 for(i=0;i<NUMVERTICES;i++){
    8095                         pe->values[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
    8096                         pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
    8097                         pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
    8098                         pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
    8099                 }
    8100         }
    8101 
    8102         /*Transform coordinate system*/
    8103         ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    8104 
    8105         /*Clean up and return*/
    8106         xDelete<int>(cs_list);
    8107         delete gauss;
    8108         return pe;
    8109 }
    8110 /*}}}*/
    8111 /*FUNCTION Penta::CreatePVectorCouplingHOFSFriction{{{*/
    8112 ElementVector* Penta::CreatePVectorCouplingHOFSFriction(void){
    8113 
    8114         /*Intermediaries*/
    8115         int         i,j;
    8116         int         approximation,analysis_type;
    8117         IssmDouble  Jdet,Jdet2d;
    8118         IssmDouble  FSreconditioning;
    8119         IssmDouble      bed_normal[3];
    8120         IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    8121         IssmDouble  viscosity, w, alpha2_gauss;
    8122         IssmDouble  dw[3];
    8123         IssmDouble      xyz_list_tria[NUMVERTICES2D][3];
    8124         IssmDouble  xyz_list[NUMVERTICES][3];
    8125         IssmDouble  basis[6]; //for the six nodes of the penta
    8126         Friction*   friction=NULL;
    8127         GaussPenta  *gauss=NULL;
    8128 
    8129         /*Initialize Element vector and return if necessary*/
    8130         if(!IsOnBed() || IsFloating()) return NULL;
    8131         inputs->GetInputValue(&approximation,ApproximationEnum);
    8132         if(approximation!=HOFSApproximationEnum) return NULL;
    8133 
    8134         int vnumnodes = this->NumberofNodesVelocity();
    8135         int pnumnodes = this->NumberofNodesPressure();
    8136         int numnodes  = vnumnodes+pnumnodes;
    8137 
    8138         /*Prepare coordinate system list*/
    8139         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    8140         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    8141         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    8142 
    8143         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,FSvelocityEnum);
    8144 
    8145         /*Retrieve all inputs and parameters*/
    8146         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8147         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    8148         this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    8149         Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
    8150         Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
    8151         Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
    8152         Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
    8153 
    8154         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    8155 
    8156         /*build friction object, used later on: */
    8157         friction=new Friction(this,3);
    8158 
    8159         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    8160         gauss=new GaussPenta(0,1,2,2);
    8161         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8162 
    8163                 gauss->GaussPoint(ig);
    8164 
    8165                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    8166                 GetNodalFunctionsP1(basis, gauss);
    8167 
    8168                 vzHO_input->GetInputValue(&w, gauss);
    8169                 vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    8170 
    8171                 NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    8172                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8173                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    8174                 friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
    8175 
    8176                 for(i=0;i<NUMVERTICES2D;i++){
    8177                         pe->values[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
    8178                         pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
    8179                         pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
    8180                 }
    8181         }
    8182 
    8183         /*Transform coordinate system*/
    8184         ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    8185 
    8186         /*Clean up and return*/
    8187         xDelete<int>(cs_list);
    8188         delete gauss;
    8189         delete friction;
    8190         return pe;
    8191 }
    8192 /*}}}*/
    8193 /*FUNCTION Penta::CreatePVectorStressbalanceHoriz{{{*/
    8194 ElementVector* Penta::CreatePVectorStressbalanceHoriz(void){
    8195 
    8196         int approximation;
    8197         inputs->GetInputValue(&approximation,ApproximationEnum);
    8198 
    8199         switch(approximation){
    8200                 case SSAApproximationEnum:
    8201                         return CreatePVectorStressbalanceSSA();
    8202                 case HOApproximationEnum:
    8203                         return CreatePVectorStressbalanceHO();
    8204                 case L1L2ApproximationEnum:
    8205                         return CreatePVectorStressbalanceL1L2();
    8206                 case SIAApproximationEnum:
    8207                         return NULL;
    8208                 case NoneApproximationEnum:
    8209                         return NULL;
    8210                 case FSApproximationEnum:
    8211                         return CreatePVectorStressbalanceFS();
    8212                 case SSAHOApproximationEnum:
    8213                         return CreatePVectorStressbalanceSSAHO();
    8214                 case SSAFSApproximationEnum:
    8215                         return CreatePVectorStressbalanceSSAFS();
    8216                 case HOFSApproximationEnum:
    8217                         return CreatePVectorStressbalanceHOFS();
    8218                 default:
    8219                         _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    8220         }
    8221 }
    8222 /*}}}*/
    8223 /*FUNCTION Penta::CreatePVectorStressbalanceSSAHO{{{*/
    8224 ElementVector* Penta::CreatePVectorStressbalanceSSAHO(void){
    8225 
    8226         /*compute all load vectors for this element*/
    8227         ElementVector* pe1=CreatePVectorStressbalanceSSA();
    8228         ElementVector* pe2=CreatePVectorStressbalanceHO();
    8229         ElementVector* pe =new ElementVector(pe1,pe2);
    8230 
    8231         /*clean-up and return*/
    8232         delete pe1;
    8233         delete pe2;
    8234         return pe;
    8235 }
    8236 /*}}}*/
    8237 /*FUNCTION Penta::CreatePVectorStressbalanceSSAFS{{{*/
    8238 ElementVector* Penta::CreatePVectorStressbalanceSSAFS(void){
    8239 
    8240         /*compute all load vectors for this element*/
    8241         int init = this->element_type;
    8242         this->element_type=P1Enum;
    8243         ElementVector* pe1=CreatePVectorStressbalanceSSA();
    8244         this->element_type=init;
    8245         ElementVector* pe2=CreatePVectorStressbalanceFS();
    8246         int indices[3]={18,19,20};
    8247         this->element_type=MINIcondensedEnum;
    8248         ElementMatrix* Ke = CreateKMatrixStressbalanceFS();
    8249         this->element_type=init;
    8250         pe2->StaticCondensation(Ke,3,&indices[0]);
    8251         delete Ke;
    8252         ElementVector* pe3=CreatePVectorCouplingSSAFS();
    8253         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    8254 
    8255         /*clean-up and return*/
    8256         delete pe1;
    8257         delete pe2;
    8258         delete pe3;
    8259         return pe;
    8260 }
    8261 /*}}}*/
    8262 /*FUNCTION Penta::CreatePVectorStressbalanceHOFS{{{*/
    8263 ElementVector* Penta::CreatePVectorStressbalanceHOFS(void){
    8264 
    8265         /*compute all load vectors for this element*/
    8266         int init = this->element_type;
    8267         this->element_type=P1Enum;
    8268         ElementVector* pe1=CreatePVectorStressbalanceHO();
    8269         this->element_type=init;
    8270         ElementVector* pe2=CreatePVectorStressbalanceFS();
    8271         int indices[3]={18,19,20};
    8272         this->element_type=MINIcondensedEnum;
    8273         ElementMatrix* Ke = CreateKMatrixStressbalanceFS();
    8274         this->element_type=init;
    8275         pe2->StaticCondensation(Ke,3,&indices[0]);
    8276         delete Ke;
    8277         ElementVector* pe3=CreatePVectorCouplingHOFS();
    8278         ElementVector* pe =new ElementVector(pe1,pe2,pe3);
    8279 
    8280         /*clean-up and return*/
    8281         delete pe1;
    8282         delete pe2;
    8283         delete pe3;
    8284         return pe;
    8285 }
    8286 /*}}}*/
    8287 /*FUNCTION Penta::CreatePVectorStressbalanceSIA{{{*/
    8288 ElementVector* Penta::CreatePVectorStressbalanceSIA(void){
    8289 
    8290         /*Intermediaries*/
    8291         int         i,j;
    8292         int         node0,node1;
    8293         IssmDouble  connectivity[2];
    8294         IssmDouble  Jdet;
    8295         IssmDouble  xyz_list[NUMVERTICES][3];
    8296         IssmDouble  xyz_list_segment[2][3];
    8297         IssmDouble  z_list[NUMVERTICES];
    8298         IssmDouble  slope[2];
    8299         IssmDouble  slope2,constant_part,drag;
    8300         IssmDouble  rho_ice,gravity,n,B;
    8301         IssmDouble  ub,vb,z_g,surface,thickness;
    8302         GaussPenta* gauss=NULL;
    8303 
    8304         /*Fetch number of nodes and dof for this finite element*/
    8305         int numnodes = this->NumberofNodes(); _assert_(numnodes==6);
    8306 
    8307         /*Initialize Element vector*/
    8308         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters);
    8309 
    8310         /*Retrieve all inputs and parameters*/
    8311         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8312         rho_ice=matpar->GetRhoIce();
    8313         gravity=matpar->GetG();
    8314         n=material->GetN();
    8315         B=material->GetB();
    8316         Input* thickness_input=inputs->GetInput(ThicknessEnum);  _assert_(thickness_input);
    8317         Input* surface_input=inputs->GetInput(SurfaceEnum);      _assert_(surface_input);
    8318         Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
    8319         Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
    8320         Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
    8321         for(i=0;i<NUMVERTICES;i++)z_list[i]=xyz_list[i][2];
    8322 
    8323         /*Loop on the three segments*/
    8324         for(i=0;i<3;i++){
    8325 
    8326                 node0=i; node1=i+3;
    8327 
    8328                 for(j=0;j<3;j++){
    8329                         xyz_list_segment[0][j]=xyz_list[node0][j];
    8330                         xyz_list_segment[1][j]=xyz_list[node1][j];
    8331                 }
    8332 
    8333                 connectivity[0]=(IssmDouble)vertices[node0]->Connectivity();
    8334                 connectivity[1]=(IssmDouble)vertices[node1]->Connectivity();
    8335 
    8336                 /*Loop on the Gauss points: */
    8337                 gauss=new GaussPenta(node0,node1,3);
    8338                 for(int ig=gauss->begin();ig<gauss->end();ig++){
    8339                         gauss->GaussPoint(ig);
    8340 
    8341                         slopex_input->GetInputValue(&slope[0],gauss);
    8342                         slopey_input->GetInputValue(&slope[1],gauss);
    8343                         surface_input->GetInputValue(&surface,gauss);
    8344                         thickness_input->GetInputValue(&thickness,gauss);
    8345 
    8346                         slope2=pow(slope[0],2)+pow(slope[1],2);
    8347                         constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
    8348 
    8349                         PentaRef::GetInputValue(&z_g,&z_list[0],gauss);
    8350                         GetSegmentJacobianDeterminant(&Jdet,&xyz_list_segment[0][0],gauss);
    8351 
    8352                         if(IsOnSurface()){
    8353                                 pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight/connectivity[1];
    8354                                 pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight/connectivity[1];
    8355                         }
    8356                         else{/*connectivity is too large, should take only half on it*/
    8357                                 pe->values[2*node1+0]+=constant_part*pow((surface-z_g)/B,n)*slope[0]*Jdet*gauss->weight*2./connectivity[1];
    8358                                 pe->values[2*node1+1]+=constant_part*pow((surface-z_g)/B,n)*slope[1]*Jdet*gauss->weight*2./connectivity[1];
    8359                         }
    8360                 }
    8361 
    8362                 /*Deal with lower surface*/
    8363                 if(IsOnBed()){
    8364                         drag_input->GetInputValue(&drag,gauss);
    8365 
    8366                         /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
    8367                         ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
    8368                         vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
    8369                         ///*Ritz et al. 1996*/
    8370                         //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
    8371                         //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
    8372                         ///*Rutt et al. 2009*/
    8373                         //ub=-drag*rho_ice*gravity*thickness*slope[0];
    8374                         //vb=-drag*rho_ice*gravity*thickness*slope[1];
    8375 
    8376                         pe->values[2*node0+0]+=ub/connectivity[0];
    8377                         pe->values[2*node0+1]+=vb/connectivity[0];
    8378                 }
    8379                 delete gauss;
    8380         }
    8381 
    8382         /*Clean up and return*/
    8383         return pe;
    8384 }
    8385 /*}}}*/
    8386 /*FUNCTION Penta::CreatePVectorStressbalanceSSA{{{*/
    8387 ElementVector* Penta::CreatePVectorStressbalanceSSA(void){
    8388 
    8389         if (!IsOnBed()) return NULL;
    8390 
    8391         /*Call Tria function*/
    8392         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    8393         ElementVector* pe=tria->CreatePVectorStressbalanceSSA();
    8394         delete tria->material; delete tria;
    8395 
    8396         /*Clean up and return*/
    8397         return pe;
    8398 }
    8399 /*}}}*/
    8400 /*FUNCTION Penta::CreatePVectorStressbalanceL1L2{{{*/
    8401 ElementVector* Penta::CreatePVectorStressbalanceL1L2(void){
    8402 
    8403         if (!IsOnBed()) return NULL;
    8404 
    8405         /*Call Tria function*/
    8406         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    8407         ElementVector* pe=tria->CreatePVectorStressbalanceSSA();
    8408         delete tria->material; delete tria;
    8409 
    8410         /*Clean up and return*/
    8411         return pe;
    8412 }
    8413 /*}}}*/
    8414 /*FUNCTION Penta::CreatePVectorStressbalanceHO{{{*/
    8415 ElementVector* Penta::CreatePVectorStressbalanceHO(void){
    8416 
    8417         /*compute all load vectors for this element*/
    8418         ElementVector* pe1=CreatePVectorStressbalanceHODrivingStress();
    8419         ElementVector* pe2=CreatePVectorStressbalanceHOFront();
    8420         ElementVector* pe =new ElementVector(pe1,pe2);
    8421 
    8422         /*clean-up and return*/
    8423         delete pe1;
    8424         delete pe2;
    8425         return pe;
    8426 }
    8427 /*}}}*/
    8428 /*FUNCTION Penta::CreatePVectorStressbalanceHODrivingStress{{{*/
    8429 ElementVector* Penta::CreatePVectorStressbalanceHODrivingStress(void){
    8430 
    8431         /*Intermediaries*/
    8432         IssmDouble  Jdet;
    8433         IssmDouble  slope[3]; //do not put 2! this goes into GetInputDerivativeValue, which addresses slope[3] also!
    8434         IssmDouble  driving_stress_baseline,thickness;
    8435         IssmDouble  xyz_list[NUMVERTICES][3];
    8436         GaussPenta  *gauss=NULL;
    8437 
    8438         /*Fetch number of nodes and dof for this finite element*/
    8439         int numnodes = this->NumberofNodes();
    8440 
    8441         /*Initialize vector*/
    8442         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
    8443         IssmDouble*    basis = xNewZeroInit<IssmDouble>(numnodes);
    8444 
    8445         /*Retrieve all inputs and parameters*/
    8446         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8447         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    8448         Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
    8449 
    8450         /* Start  looping on the number of gaussian points: */
    8451         gauss=new GaussPenta(3,3);
    8452         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8453 
    8454                 gauss->GaussPoint(ig);
    8455 
    8456                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8457                 GetNodalFunctions(basis,gauss);
    8458 
    8459                 thickness_input->GetInputValue(&thickness, gauss);
    8460                 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    8461 
    8462                 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
    8463 
    8464                 for(int i=0;i<numnodes;i++){
    8465                         pe->values[i*NDOF2+0]+= -driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
    8466                         pe->values[i*NDOF2+1]+= -driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
    8467                 }
    8468         }
    8469 
    8470         /*Transform coordinate system*/
    8471         ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
    8472 
    8473         /*Clean up and return*/
    8474         xDelete<IssmDouble>(basis);
    8475         delete gauss;
    8476         return pe;
    8477 }
    8478 /*}}}*/
    8479 /*FUNCTION Penta::CreatePVectorStressbalanceHOFront{{{*/
    8480 ElementVector* Penta::CreatePVectorStressbalanceHOFront(void){
    8481 
    8482         /*If no front, return NULL*/
    8483         if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
    8484 
    8485         /*Intermediaries*/
    8486         IssmDouble  rho_ice,rho_water,gravity;
    8487         IssmDouble  Jdet,surface,z_g,water_pressure,ice_pressure,air_pressure;
    8488         IssmDouble  surface_under_water,base_under_water,pressure;
    8489         IssmDouble  xyz_list[NUMVERTICES][3];
    8490         IssmDouble* xyz_list_front = NULL;
    8491         IssmDouble  area_coordinates[4][3];
    8492         IssmDouble  normal[3];
    8493         GaussPenta* gauss;
    8494 
    8495         /*Fetch number of nodes and dof for this finite element*/
    8496         int numnodes = this->NumberofNodes();
    8497 
    8498         /*Initialize Element vector and other vectors*/
    8499         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8500         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
    8501         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    8502 
    8503         /*Retrieve all inputs and parameters*/
    8504         Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
    8505         rho_water=matpar->GetRhoWater();
    8506         rho_ice  =matpar->GetRhoIce();
    8507         gravity  =matpar->GetG();
    8508         ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    8509         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
    8510         NormalSection(&normal[0],xyz_list_front);
    8511 
    8512         /*Initialize gauss points*/
    8513         IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
    8514         IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
    8515         if(zmax>0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,10); //refined in vertical because of the sea level discontinuity
    8516         else                 gauss=new GaussPenta(area_coordinates,3,3);
    8517 
    8518         /* Start  looping on the number of gaussian points: */
    8519         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8520 
    8521                 gauss->GaussPoint(ig);
    8522                 surface_input->GetInputValue(&surface,gauss);
    8523                 z_g=GetZcoord(gauss);
    8524                 GetNodalFunctions(basis,gauss);
    8525                 GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss);
    8526 
    8527                 water_pressure = rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
    8528                 ice_pressure   = rho_ice*gravity*(surface-z_g);
    8529                 air_pressure   = 0;
    8530                 pressure       = ice_pressure + water_pressure + air_pressure;
    8531 
    8532                 for (int i=0;i<numnodes;i++){
    8533                         pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
    8534                         pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
    8535                 }
    8536         }
    8537 
    8538         /*Transform coordinate system*/
    8539         ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
    8540 
    8541         /*Clean up and return*/
    8542         xDelete<IssmDouble>(basis);
    8543         xDelete<IssmDouble>(xyz_list_front);
    8544         delete gauss;
    8545         return pe;
    8546 }
    8547 /*}}}*/
    8548 /*FUNCTION Penta::CreatePVectorStressbalanceFS {{{*/
    8549 ElementVector* Penta::CreatePVectorStressbalanceFS(void){
    8550 
    8551         ElementVector* pe1;
    8552         ElementVector* pe2;
    8553         ElementVector* pe3;
    8554         ElementVector* pe;
    8555 
    8556         /*compute all stiffness matrices for this element*/
    8557         pe1=CreatePVectorStressbalanceFSViscous();
    8558         pe2=CreatePVectorStressbalanceFSShelf();
    8559         pe3=CreatePVectorStressbalanceFSFront();
    8560         pe =new ElementVector(pe1,pe2,pe3);
    8561 
    8562         /*clean-up and return*/
    8563         delete pe1;
    8564         delete pe2;
    8565         delete pe3;
    8566         return pe;
    8567 }
    8568 /*}}}*/
    8569 /*FUNCTION Penta::CreatePVectorStressbalanceFSFront{{{*/
    8570 ElementVector* Penta::CreatePVectorStressbalanceFSFront(void){
    8571 
    8572         /*If no front, return NULL*/
    8573         if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
    8574 
    8575         /*Fetch number of nodes and dof for this finite element*/
    8576         int         i;
    8577         IssmDouble  rho_ice,rho_water,gravity;
    8578         IssmDouble  Jdet,z_g,water_pressure,air_pressure;
    8579         IssmDouble  surface_under_water,base_under_water,pressure;
    8580         IssmDouble  xyz_list[NUMVERTICES][3];
    8581         IssmDouble* xyz_list_front = NULL;
    8582         IssmDouble  area_coordinates[4][3];
    8583         IssmDouble  normal[3];
    8584         GaussPenta* gauss;
    8585 
    8586         /*Fetch number of nodes and dof for this finite element*/
    8587         int vnumnodes = this->NumberofNodesVelocity();
    8588         int pnumnodes = this->NumberofNodesPressure();
    8589 
    8590         /*Prepare coordinate system list*/
    8591         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    8592         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    8593         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    8594 
    8595         /*Initialize Element matrix and vectors*/
    8596         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8597         ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    8598         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    8599 
    8600         /*Retrieve all inputs and parameters*/
    8601         rho_water=matpar->GetRhoWater();
    8602         rho_ice  =matpar->GetRhoIce();
    8603         gravity  =matpar->GetG();
    8604         ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    8605         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],4);
    8606         NormalSection(&normal[0],xyz_list_front);
    8607 
    8608         /*Start looping on Gaussian points*/
    8609         IssmDouble zmax=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
    8610         IssmDouble zmin=xyz_list[0][2]; for(int i=1;i<6;i++) if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
    8611         if(zmax>0 && zmin<0) gauss=new GaussPenta(area_coordinates,3,30); //refined in vertical because of the sea level discontinuity
    8612         else                 gauss=new GaussPenta(area_coordinates,3,3);
    8613 
    8614         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8615 
    8616                 gauss->GaussPoint(ig);
    8617                 z_g=GetZcoord(gauss);
    8618                 GetNodalFunctionsVelocity(vbasis,gauss);
    8619                 GetQuadJacobianDeterminant(&Jdet,xyz_list_front,gauss);
    8620 
    8621                 water_pressure=rho_water*gravity*min(0.,z_g);//0 if the gaussian point is above water level
    8622                 air_pressure=0;
    8623                 pressure = water_pressure + air_pressure;
    8624 
    8625                 for(i=0;i<vnumnodes;i++){
    8626                         pe->values[3*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
    8627                         pe->values[3*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
    8628                         pe->values[3*i+2]+= pressure*Jdet*gauss->weight*normal[2]*vbasis[i];
    8629                 }
    8630         }
    8631 
    8632         /*Transform coordinate system*/
    8633         ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    8634 
    8635         /*Clean up and return*/
    8636         xDelete<int>(cs_list);
    8637         xDelete<IssmDouble>(xyz_list_front);
    8638         xDelete<IssmDouble>(vbasis);
    8639         delete gauss;
    8640         return pe;
    8641 }
    8642 /*}}}*/
    8643 /*FUNCTION Penta::PVectorGLSstabilization{{{*/
    8644 void Penta::PVectorGLSstabilization(ElementVector* pe){
    8645 
    8646         /*Intermediaries*/
    8647         int        i,j;
    8648         IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab,viscosity;
    8649         IssmDouble forcex,forcey,forcez,diameter,FSreconditioning;
    8650         IssmDouble xyz_list[NUMVERTICES][3];
    8651         IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    8652         GaussPenta *gauss=NULL;
    8653 
    8654         /*Stabilization*/
    8655         IssmDouble dbasis[3][6];
    8656         IssmDouble dmu[3];
    8657         IssmDouble mu;
    8658         IssmDouble dnodalbasis[6][6][3];
    8659         IssmDouble SW[6][4][4];
    8660         int p,ii;
    8661         int c=3; //index of pressure
    8662 
    8663         /*Retrieve all inputs and parameters*/
    8664         rho_ice=matpar->GetRhoIce();
    8665         gravity=matpar->GetG();
    8666         B=material->GetB();
    8667         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8668         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    8669         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    8670         Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input);
    8671 
    8672         /*Find minimal length*/
    8673         diameter=MinEdgeLength(&xyz_list[0][0]);
    8674 
    8675                 gauss=new GaussPenta();
    8676                 for(int iv=0;iv<6;iv++){
    8677                         gauss->GaussVertex(iv);
    8678                         GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    8679                         for(i=0;i<6;i++){
    8680                                 for(j=0;j<3;j++){
    8681                                         dnodalbasis[i][iv][j] = dbasis[j][i];
    8682                                 }
    8683                         }
    8684                 }
    8685                 delete gauss;
    8686 
    8687         /* Start  looping on the number of gaussian points: */
    8688         gauss=new GaussPenta(5,5);
    8689         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8690 
    8691                 gauss->GaussPoint(ig);
    8692 
    8693                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8694                 this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    8695                 material->GetViscosity3dFS(&viscosity,&epsilon[0]);
    8696 
    8697                 GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
    8698                 dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
    8699                 mu = 2.*viscosity;
    8700                 for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
    8701                         SW[p][i][j]=0.;
    8702                 }
    8703                 for(p=0;p<6;p++){
    8704                         for(i=0;i<3;i++){
    8705                                 SW[p][c][i] += dbasis[i][p];
    8706                                 for(j=0;j<3;j++){
    8707                                         SW[p][i][i] += -dmu[j]*dbasis[j][p];
    8708                                         SW[p][j][i] += -dmu[j]*dbasis[i][p];
    8709                                         for(ii=0;ii<6;ii++){
    8710                                                 SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
    8711                                                 SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
    8712                                         }
    8713                                 }
    8714                         }
    8715                 }
    8716                 IssmDouble tau = 1./3.*pow(diameter,2)/(8.*2.*viscosity);
    8717                 for(p=0;p<6;p++){
    8718                         for(j=0;j<4;j++){
    8719                                 pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcex*SW[p][j][0];
    8720                                 pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*forcey*SW[p][j][1];
    8721                                 pe->values[p*4+j] += gauss->weight*Jdet*tau*rho_ice*(forcez-gravity)*SW[p][j][2];
    8722                         }
    8723                 }
    8724         }
    8725 
    8726         /*Clean up*/
    8727         delete gauss;
    8728 }
    8729 /*}}}*/
    8730 /*FUNCTION Penta::CreatePVectorStressbalanceFSViscous {{{*/
    8731 ElementVector* Penta::CreatePVectorStressbalanceFSViscous(void){
    8732 
    8733         /*Intermediaries*/
    8734         int        i;
    8735         int        approximation;
    8736         IssmDouble Jdet,gravity,rho_ice;
    8737         IssmDouble forcex,forcey,forcez;
    8738         IssmDouble xyz_list[NUMVERTICES][3];
    8739         GaussPenta *gauss=NULL;
    8740 
    8741         /*Initialize Element vector and return if necessary*/
    8742         inputs->GetInputValue(&approximation,ApproximationEnum);
    8743         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    8744 
    8745         /*Fetch number of nodes and dof for this finite element*/
    8746         int vnumnodes = this->NumberofNodesVelocity();
    8747         int pnumnodes = this->NumberofNodesPressure();
    8748 
    8749         /*Prepare coordinate system list*/
    8750         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    8751         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    8752         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    8753 
    8754         /*Initialize Element matrix and vectors*/
    8755         ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    8756         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    8757 
    8758         /*Retrieve all inputs and parameters*/
    8759         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8760         Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
    8761         Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
    8762         Input* loadingforcez_input=inputs->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
    8763         rho_ice = matpar->GetRhoIce();
    8764         gravity = matpar->GetG();
    8765 
    8766         /* Start  looping on the number of gaussian points: */
    8767         gauss=new GaussPenta(5,5);
    8768         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8769 
    8770                 gauss->GaussPoint(ig);
    8771 
    8772                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8773                 GetNodalFunctionsVelocity(vbasis, gauss);
    8774 
    8775                 loadingforcex_input->GetInputValue(&forcex,gauss);
    8776                 loadingforcey_input->GetInputValue(&forcey,gauss);
    8777                 loadingforcez_input->GetInputValue(&forcez,gauss);
    8778 
    8779                 for(i=0;i<vnumnodes;i++){
    8780                         pe->values[i*NDOF3+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
    8781                         pe->values[i*NDOF3+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
    8782                         pe->values[i*NDOF3+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
    8783                         pe->values[i*NDOF3+2] += +rho_ice*forcez *Jdet*gauss->weight*vbasis[i];
    8784                 }
    8785         }
    8786 
    8787         /*Transform coordinate system*/
    8788         ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    8789 
    8790         /*Clean up and return*/
    8791         delete gauss;
    8792         xDelete<int>(cs_list);
    8793         xDelete<IssmDouble>(vbasis);
    8794         return pe;
    8795 }
    8796 /*}}}*/
    8797 /*FUNCTION Penta::CreatePVectorStressbalanceFSShelf{{{*/
    8798 ElementVector* Penta::CreatePVectorStressbalanceFSShelf(void){
    8799 
    8800         /*Intermediaries*/
    8801         int         i,j;
    8802         int         approximation,shelf_dampening;
    8803         IssmDouble  gravity,rho_water,bed,water_pressure;
    8804         IssmDouble  damper,normal_vel,vx,vy,vz,dt;
    8805         IssmDouble      xyz_list_tria[NUMVERTICES2D][3];
    8806         IssmDouble  xyz_list[NUMVERTICES][3];
    8807         IssmDouble      bed_normal[3];
    8808         IssmDouble  dz[3];
    8809         IssmDouble  Jdet2d;
    8810         GaussPenta  *gauss=NULL;
    8811 
    8812         /*Initialize Element vector and return if necessary*/
    8813         if(!IsOnBed() || !IsFloating()) return NULL;
    8814         inputs->GetInputValue(&approximation,ApproximationEnum);
    8815         if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
    8816 
    8817         /*Fetch number of nodes and dof for this finite element*/
    8818         int vnumnodes = this->NumberofNodesVelocity();
    8819         int pnumnodes = this->NumberofNodesPressure();
    8820 
    8821         /*Prepare coordinate system list*/
    8822         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    8823         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
    8824         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    8825 
    8826         /*Initialize Element matrix and vectors*/
    8827         ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    8828         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    8829 
    8830         /*Retrieve all inputs and parameters*/
    8831         this->parameters->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
    8832         rho_water=matpar->GetRhoWater();
    8833         gravity=matpar->GetG();
    8834         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8835         Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input);
    8836         Input* vx_input=inputs->GetInput(VxEnum);   _assert_(vx_input);
    8837         Input* vy_input=inputs->GetInput(VyEnum);   _assert_(vy_input);
    8838         Input* vz_input=inputs->GetInput(VzEnum);   _assert_(vz_input);
    8839 
    8840         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    8841 
    8842         /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    8843         gauss=new GaussPenta(0,1,2,5);
    8844         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8845 
    8846                 gauss->GaussPoint(ig);
    8847 
    8848                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    8849                 GetNodalFunctionsVelocity(vbasis, gauss);
    8850 
    8851                 NormalBase(&bed_normal[0],&xyz_list_tria[0][0]);
    8852                 bed_input->GetInputValue(&bed, gauss);
    8853                 if(shelf_dampening){ //add dampening to avoid too high vertical velocities when not in hydrostatic equilibrium
    8854                         bed_input->GetInputDerivativeValue(&dz[0],&xyz_list[0][0],gauss);
    8855                         vx_input->GetInputValue(&vx, gauss);
    8856                         vy_input->GetInputValue(&vy, gauss);
    8857                         vz_input->GetInputValue(&vz, gauss);
    8858                         dt=0.;
    8859                         normal_vel=bed_normal[0]*vx+bed_normal[1]*vy+bed_normal[2]*vz;
    8860                         damper=gravity*rho_water*pow(1+pow(dz[0],2)+pow(dz[1],2),0.5)*normal_vel*dt;
    8861                 }
    8862                 else damper=0.;
    8863                 water_pressure=gravity*rho_water*bed;
    8864 
    8865                 for(i=0;i<vnumnodes;i++){
    8866                         pe->values[3*i+0]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[0];
    8867                         pe->values[3*i+1]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[1];
    8868                         pe->values[3*i+2]+=(water_pressure+damper)*gauss->weight*Jdet2d*vbasis[i]*bed_normal[2];
    8869                 }
    8870         }
    8871 
    8872         /*Transform coordinate system*/
    8873         ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    8874 
    8875         /*Clean up and return*/
    8876         delete gauss;
    8877         xDelete<int>(cs_list);
    8878         xDelete<IssmDouble>(vbasis);
    8879         return pe;
    8880 }
    8881 /*}}}*/
    8882 /*FUNCTION Penta::CreatePVectorStressbalanceVert {{{*/
    8883 ElementVector* Penta::CreatePVectorStressbalanceVert(void){
    8884 
    8885         /*compute all load vectors for this element*/
    8886         ElementVector* pe1=CreatePVectorStressbalanceVertVolume();
    8887         ElementVector* pe2=CreatePVectorStressbalanceVertBase();
    8888         ElementVector* pe =new ElementVector(pe1,pe2);
    8889 
    8890         /*clean-up and return*/
    8891         delete pe1;
    8892         delete pe2;
    8893         return pe;
    8894 }
    8895 /*}}}*/
    8896 /*FUNCTION Penta::CreatePVectorStressbalanceVertVolume {{{*/
    8897 ElementVector* Penta::CreatePVectorStressbalanceVertVolume(void){
    8898 
    8899         /*Constants*/
    8900         const int  numdof=NDOF1*NUMVERTICES;
    8901 
    8902         /*Intermediaries*/
    8903         int        approximation;
    8904         IssmDouble     Jdet;
    8905         IssmDouble     xyz_list[NUMVERTICES][3];
    8906         IssmDouble     dudx,dvdy,dwdz;
    8907         IssmDouble     du[3],dv[3],dw[3];
    8908         IssmDouble     basis[6];
    8909         GaussPenta *gauss=NULL;
    8910 
    8911         /*Initialize Element vector*/
    8912         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    8913 
    8914         /*Retrieve all inputs and parameters*/
    8915         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8916         inputs->GetInputValue(&approximation,ApproximationEnum);
    8917         Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
    8918         Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
    8919         Input* vzFS_input=NULL;
    8920         if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
    8921                 vzFS_input=inputs->GetInput(VzFSEnum); _assert_(vzFS_input);
    8922         }
    8923 
    8924         /* Start  looping on the number of gaussian points: */
    8925         gauss=new GaussPenta(2,2);
    8926         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8927 
    8928                 gauss->GaussPoint(ig);
    8929 
    8930                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8931                 GetNodalFunctionsP1(basis, gauss);
    8932 
    8933                 vx_input->GetInputDerivativeValue(&du[0],&xyz_list[0][0],gauss);
    8934                 vy_input->GetInputDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
    8935                 if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
    8936                         vzFS_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
    8937                         dwdz=dw[2];
    8938                 }
    8939                 else dwdz=0;
    8940                 dudx=du[0];
    8941                 dvdy=dv[1];
    8942 
    8943                 for(int i=0;i<numdof;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i];
    8944         }
    8945 
    8946         /*Clean up and return*/
    8947         delete gauss;
    8948         return pe;
    8949 }
    8950 /*}}}*/
    8951 /*FUNCTION Penta::CreatePVectorStressbalanceVertBase {{{*/
    8952 ElementVector* Penta::CreatePVectorStressbalanceVertBase(void){
    8953 
    8954         /*Constants*/
    8955         const int    numdof=NDOF1*NUMVERTICES;
    8956 
    8957         /*Intermediaries */
    8958         int        i,j;
    8959         int        approximation;
    8960         IssmDouble xyz_list[NUMVERTICES][3];
    8961         IssmDouble xyz_list_tria[NUMVERTICES2D][3];
    8962         IssmDouble Jdet2d;
    8963         IssmDouble vx,vy,vz,dbdx,dbdy,basalmeltingvalue;
    8964         IssmDouble slope[3];
    8965         IssmDouble basis[NUMVERTICES];
    8966         GaussPenta* gauss=NULL;
    8967 
    8968         if (!IsOnBed()) return NULL;
    8969 
    8970         /*Initialize Element vector*/
    8971         ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
    8972 
    8973         /*Retrieve all inputs and parameters*/
    8974         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8975         for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
    8976         inputs->GetInputValue(&approximation,ApproximationEnum);
    8977         Input* bed_input=inputs->GetInput(BedEnum);                                _assert_(bed_input);
    8978         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
    8979         Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    8980         Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    8981         Input* vzFS_input=NULL;
    8982         if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
    8983                 vzFS_input=inputs->GetInput(VzFSEnum);       _assert_(vzFS_input);
    8984         }
    8985 
    8986         /* Start  looping on the number of gaussian points: */
    8987         gauss=new GaussPenta(0,1,2,2);
    8988         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8989 
    8990                 gauss->GaussPoint(ig);
    8991 
    8992                 basal_melting_input->GetInputValue(&basalmeltingvalue, gauss);
    8993                 bed_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    8994                 vx_input->GetInputValue(&vx, gauss);
    8995                 vy_input->GetInputValue(&vy, gauss);
    8996                 if(approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
    8997                         vzFS_input->GetInputValue(&vz, gauss);
    8998                 }
    8999                 else vz=0;
    9000 
    9001                 dbdx=slope[0];
    9002                 dbdy=slope[1];
    9003 
    9004                 GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
    9005                 GetNodalFunctionsP1(&basis[0], gauss);
    9006 
    9007                 for(i=0;i<numdof;i++) pe->values[i]+=-Jdet2d*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
    9008         }
    9009 
    9010         /*Clean up and return*/
    9011         delete gauss;
    9012         return pe;
    9013 }
    9014 /*}}}*/
    90156803/*FUNCTION Penta::CreateJacobianStressbalanceHoriz{{{*/
    90166804ElementMatrix* Penta::CreateJacobianStressbalanceHoriz(void){
    … …  
    92467034}
    92477035/*}}}*/
    9248 /*FUNCTION Penta::CreatePVectorBalancethickness {{{*/
    9249 ElementVector* Penta::CreatePVectorBalancethickness(void){
    9250 
    9251         if (!IsOnBed()) return NULL;
    9252 
    9253         /*Depth Averaging Vx and Vy*/
    9254         this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
    9255         this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
    9256 
    9257         /*Call Tria function*/
    9258         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    9259         ElementVector* pe=tria->CreatePVectorBalancethickness();
    9260         delete tria->material; delete tria;
    9261 
    9262         /*Delete Vx and Vy averaged*/
    9263         this->inputs->DeleteInput(VxAverageEnum);
    9264         this->inputs->DeleteInput(VyAverageEnum);
    9265 
    9266         /*Clean up and return*/
    9267         return pe;
    9268 }
    9269 /*}}}*/
    92707036#endif
    92717037
    … …  
    93087074        /*clean up and return*/
    93097075        return Ke;
    9310 }
    9311 /*}}}*/
    9312 /*FUNCTION Penta::CreatePVectorHydrologyDCInefficient {{{*/
    9313 ElementVector* Penta::CreatePVectorHydrologyDCInefficient(void){
    9314 
    9315         if (!IsOnBed()) return NULL;
    9316 
    9317         /*Call Tria function*/
    9318         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    9319         ElementVector* pe=tria->CreatePVectorHydrologyDCInefficient();
    9320         delete tria->material; delete tria;
    9321 
    9322         /*Clean up and return*/
    9323         return pe;
    9324 }
    9325 /*}}}*/
    9326 /*FUNCTION Penta::CreatePVectorHydrologyDCEfficient {{{*/
    9327 ElementVector* Penta::CreatePVectorHydrologyDCEfficient(void){
    9328 
    9329         if (!IsOnBed()) return NULL;
    9330 
    9331         /*Call Tria function*/
    9332         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    9333         ElementVector* pe=tria->CreatePVectorHydrologyDCEfficient();
    9334         delete tria->material; delete tria;
    9335 
    9336         /*Clean up and return*/
    9337         return pe;
    9338 }
    9339 /*}}}*/
    9340 /*FUNCTION Penta::CreatePVectorL@ProjectionEPL {{{*/
    9341 ElementVector* Penta::CreatePVectorL2ProjectionEPL(void){
    9342 
    9343         if (!IsOnBed()) return NULL;
    9344 
    9345         /*Call Tria function*/
    9346         Tria* tria=(Tria*)SpawnTria(0); //lower face is 0, upper face is 1.
    9347         ElementVector* pe=tria->CreatePVectorL2ProjectionEPL();
    9348         delete tria->material; delete tria;
    9349 
    9350         /*Clean up and return*/
    9351         return pe;
    93527076}
    93537077/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Penta.h

    r16956 r16982  
    199199                ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
    200200                ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    201                 ElementVector* CreatePVector(void);
    202                 ElementVector* CreatePVectorMasstransport(void);
    203                 ElementVector* CreatePVectorFreeSurfaceTop(void);
    204                 ElementVector* CreatePVectorFreeSurfaceBase(void);
    205                 ElementVector* CreatePVectorL2ProjectionBase(void);
    206201                IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);
    207202                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);
    … …  
    301296                ElementMatrix* CreateJacobianStressbalanceHO(void);
    302297                ElementMatrix* CreateJacobianStressbalanceFS(void);
    303                 ElementVector* CreatePVectorCouplingSSAFS(void);
    304                 ElementVector* CreatePVectorCouplingSSAFSViscous(void);
    305                 ElementVector* CreatePVectorCouplingSSAFSFriction(void);
    306                 ElementVector* CreatePVectorCouplingHOFS(void);
    307                 ElementVector* CreatePVectorCouplingHOFSViscous(void);
    308                 ElementVector* CreatePVectorCouplingHOFSFriction(void);
    309                 ElementVector* CreatePVectorStressbalanceHoriz(void);
    310                 ElementVector* CreatePVectorStressbalanceSIA(void);
    311                 ElementVector* CreatePVectorStressbalanceSSA(void);
    312                 ElementVector* CreatePVectorStressbalanceSSAHO(void);
    313                 ElementVector* CreatePVectorStressbalanceSSAFS(void);
    314                 ElementVector* CreatePVectorStressbalanceL1L2(void);
    315                 ElementVector* CreatePVectorStressbalanceHO(void);
    316                 ElementVector* CreatePVectorStressbalanceHODrivingStress(void);
    317                 ElementVector* CreatePVectorStressbalanceHOFront(void);
    318                 ElementVector* CreatePVectorStressbalanceHOFS(void);
    319                 ElementVector* CreatePVectorStressbalanceFS(void);
    320                 ElementVector* CreatePVectorStressbalanceFSFront(void);
    321                 ElementVector* CreatePVectorStressbalanceFSViscous(void);
    322                 void           PVectorGLSstabilization(ElementVector* pe);
    323                 ElementVector* CreatePVectorStressbalanceFSShelf(void);
    324                 ElementVector* CreatePVectorStressbalanceVert(void);
    325                 ElementVector* CreatePVectorStressbalanceVertVolume(void);
    326                 ElementVector* CreatePVectorStressbalanceVertBase(void);
    327298                #endif
    328299
    329300                #ifdef _HAVE_CONTROL_
    330                 ElementVector* CreatePVectorAdjointHoriz(void);
    331301                ElementMatrix* CreateKMatrixAdjointSSA2d(void);
    332302                ElementMatrix* CreateKMatrixAdjointHO(void);
    333303                ElementMatrix* CreateKMatrixAdjointFS(void);
    334                 ElementVector* CreatePVectorAdjointSSA(void);
    335                 ElementVector* CreatePVectorAdjointHO(void);
    336                 ElementVector* CreatePVectorAdjointFS(void);
    337304                #endif
    338305
    … …  
    341308                ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
    342309                ElementMatrix* CreateEPLDomainMassMatrix(void);
    343                 ElementVector* CreatePVectorHydrologyDCInefficient(void);
    344                 ElementVector* CreatePVectorHydrologyDCEfficient(void);
    345                 ElementVector* CreatePVectorL2ProjectionEPL(void);
    346310                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode);
    347311                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index){_error_("not implemented yet");};
    … …  
    362326                ElementMatrix* CreateKMatrixThermalVolume(void);
    363327                ElementMatrix* CreateKMatrixThermalShelf(void);
    364                 ElementVector* CreatePVectorEnthalpy(void);
    365                 ElementVector* CreatePVectorEnthalpyVolume(void);
    366                 ElementVector* CreatePVectorEnthalpyShelf(void);
    367                 ElementVector* CreatePVectorEnthalpySheet(void);
    368                 ElementVector* CreatePVectorMelting(void);
    369                 ElementVector* CreatePVectorThermal(void);
    370                 ElementVector* CreatePVectorThermalVolume(void);
    371                 ElementVector* CreatePVectorThermalShelf(void);
    372                 ElementVector* CreatePVectorThermalSheet(void);
    373328                void           UpdateBasalConstraintsEnthalpy(void);
    374329                void           ComputeBasalMeltingrate(void);
    … …  
    378333                #ifdef _HAVE_BALANCED_
    379334                ElementMatrix* CreateKMatrixBalancethickness(void);
    380                 ElementVector* CreatePVectorBalancethickness(void);
    381335                #endif
    382336                /*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Seg.cpp

    r16913 r16982  
    333333}
    334334/*}}}*/
    335 /*FUNCTION Seg::CreatePVectorL2Projection {{{*/
    336 ElementVector* Seg::CreatePVectorL2Projection(void){
    337 
    338         /*Intermediaries */
    339         int        i,input_enum;
    340         IssmDouble Jdet;
    341         IssmDouble xyz_list[NUMVERTICES][3];
    342         IssmDouble value;
    343         Input*     input  = NULL;
    344         Input*     input2 = NULL;
    345 
    346         /*Fetch number of nodes and dof for this finite element*/
    347         int numnodes = this->NumberofNodes();
    348 
    349         /*Initialize Element vector*/
    350         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    351         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    352 
    353         /*Retrieve all inputs and parameters*/
    354         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    355         this->parameters->FindParam(&input_enum,InputToL2ProjectEnum);
    356         switch(input_enum){
    357                 case SurfaceSlopeXEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break;
    358                 case BedSlopeXEnum:     input2 = inputs->GetInput(BedEnum);     _assert_(input2); break;
    359                 default: input = inputs->GetInput(input_enum);
    360         }
    361 
    362         /* Start  looping on the number of gaussian points: */
    363         GaussSeg* gauss=new GaussSeg(2);
    364         for(int ig=gauss->begin();ig<gauss->end();ig++){
    365 
    366                 gauss->GaussPoint(ig);
    367 
    368                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    369                 GetNodalFunctions(basis,gauss);
    370 
    371                 if(input2){
    372                         input2->GetInputDerivativeValue(&value,&xyz_list[0][0],gauss);
    373                 }
    374                 else{
    375                         input->GetInputValue(&value,gauss);
    376                 }
    377 
    378                 for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
    379         }
    380 
    381         /*Clean up and return*/
    382         xDelete<IssmDouble>(basis);
    383         delete gauss;
    384         return pe;
    385 }
    386 /*}}}*/
    387 /*FUNCTION Seg::CreatePVectorFreeSurfaceTop {{{*/
    388 ElementVector* Seg::CreatePVectorFreeSurfaceTop(void){
    389 
    390         /*Intermediaries */
    391         IssmDouble Jdet,dt;
    392         IssmDouble ms,surface,vy;
    393         IssmDouble xyz_list[NUMVERTICES][3];
    394 
    395         /*Fetch number of nodes and dof for this finite element*/
    396         int numnodes = this->NumberofNodes();
    397 
    398         /*Initialize Element vector and other vectors*/
    399         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    400         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    401 
    402         /*Retrieve all inputs and parameters*/
    403         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    404         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    405         Input* vy_input     = inputs->GetInput(VyEnum);                         _assert_(vy_input);
    406         Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
    407         Input* surface_input= inputs->GetInput(SurfaceEnum);                    _assert_(surface_input);
    408 
    409         /*Initialize mb_correction to 0, do not forget!:*/
    410         /* Start  looping on the number of gaussian points: */
    411         GaussSeg* gauss=new GaussSeg(2);
    412         for(int ig=gauss->begin();ig<gauss->end();ig++){
    413 
    414                 gauss->GaussPoint(ig);
    415 
    416                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    417                 GetNodalFunctions(basis,gauss);
    418 
    419                 vy_input->GetInputValue(&vy,gauss);
    420                 ms_input->GetInputValue(&ms,gauss);
    421                 surface_input->GetInputValue(&surface,gauss);
    422 
    423                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vy)*basis[i];
    424         }
    425 
    426         /*Clean up and return*/
    427         xDelete<IssmDouble>(basis);
    428         delete gauss;
    429         return pe;
    430 }
    431 /*}}}*/
    432 /*FUNCTION Seg::CreatePVectorFreeSurfaceBase {{{*/
    433 ElementVector* Seg::CreatePVectorFreeSurfaceBase(void){
    434 
    435         /*Intermediaries */
    436         IssmDouble Jdet,dt;
    437         IssmDouble mb,mb_correction,bed,vy;
    438         IssmDouble xyz_list[NUMVERTICES][3];
    439 
    440         /*Fetch number of nodes and dof for this finite element*/
    441         int numnodes = this->NumberofNodes();
    442 
    443         /*Initialize Element vector and other vectors*/
    444         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    445         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    446 
    447         /*Retrieve all inputs and parameters*/
    448         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    449         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    450         Input* vy_input            = inputs->GetInput(VyEnum);                         _assert_(vy_input);
    451         Input* mb_input            = inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
    452         Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
    453         Input* bed_input           = inputs->GetInput(BedEnum);                        _assert_(bed_input);
    454 
    455         /*Initialize mb_correction to 0, do not forget!:*/
    456         /* Start  looping on the number of gaussian points: */
    457         GaussSeg* gauss=new GaussSeg(2);
    458         for(int ig=gauss->begin();ig<gauss->end();ig++){
    459 
    460                 gauss->GaussPoint(ig);
    461 
    462                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    463                 GetNodalFunctions(basis,gauss);
    464 
    465                 vy_input->GetInputValue(&vy,gauss);
    466                 mb_input->GetInputValue(&mb,gauss);
    467                 bed_input->GetInputValue(&bed,gauss);
    468                 if(mb_correction_input)
    469                  mb_correction_input->GetInputValue(&mb_correction,gauss);
    470                 else
    471                  mb_correction=0.;
    472 
    473                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vy)*basis[i];
    474         }
    475 
    476         /*Clean up and return*/
    477         xDelete<IssmDouble>(basis);
    478         delete gauss;
    479         return pe;
    480 }
    481 /*}}}*/
    482335/*FUNCTION Seg::GetNumberOfNodes;{{{*/
    483336int Seg::GetNumberOfNodes(void){

  • issm/trunk-jpl/src/c/classes/Elements/Seg.h

    r16956 r16982  
    7171                ElementMatrix* CreateKMatrix(void){_error_("not implemented yet");};
    7272                void        CreateDVector(Vector<IssmDouble>* df){_error_("not implemented yet");};
    73                 ElementVector* CreatePVector(void){_error_("not implemented yet");};
    74                 ElementVector* CreatePVectorL2Projection(void);
    7573                void        CreateJacobianMatrix(Matrix<IssmDouble>* Jff){_error_("not implemented yet");};
    7674                void        Delta18oParameterization(void){_error_("not implemented yet");};
    … …  
    262260                ElementMatrix* CreateKMatrixFreeSurfaceTop(void);
    263261                ElementMatrix* CreateKMatrixFreeSurfaceBase(void);
    264                 ElementVector* CreatePVectorFreeSurfaceTop(void);
    265                 ElementVector* CreatePVectorFreeSurfaceBase(void);
    266262                IssmDouble     GetSize(void);
    267263};

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r16974 r16982  
    374374}
    375375/*}}}*/
    376 /*FUNCTION Tria::CreatePVector(void){{{*/
    377 ElementVector* Tria::CreatePVector(void){
    378 
    379         /*retrive parameters: */
    380         int meshtype;
    381         int analysis_type;
    382         parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    383 
    384         /*asserts: {{{*/
    385         /*if debugging mode, check that all pointers exist*/
    386         _assert_(this->nodes && this->material && this->matpar && this->parameters && this->inputs);
    387         /*}}}*/
    388 
    389         /*Skip if water element*/
    390         if(NoIceInElement()) return NULL;
    391 
    392         /*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
    393         switch(analysis_type){
    394                 #ifdef _HAVE_STRESSBALANCE_
    395                 case StressbalanceAnalysisEnum:
    396                         int approximation;
    397                         inputs->GetInputValue(&approximation,ApproximationEnum);
    398                         switch(approximation){
    399                                 case SSAApproximationEnum:
    400                                         return CreatePVectorStressbalanceSSA();
    401                                 case FSApproximationEnum:
    402                                         return CreatePVectorStressbalanceFS();
    403                                 case SIAApproximationEnum:
    404                                         return NULL;
    405                                 case NoneApproximationEnum:
    406                                         return NULL;
    407                                 default:
    408                                         _error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
    409                         }
    410                         break;
    411                 case StressbalanceSIAAnalysisEnum:
    412                         return CreatePVectorStressbalanceSIA();
    413                         break;
    414                 #endif
    415                 #ifdef _HAVE_DAMAGE_
    416                 case DamageEvolutionAnalysisEnum:
    417                         return CreatePVectorDamageEvolution();
    418                         break;
    419                 #endif
    420                 case L2ProjectionBaseAnalysisEnum:
    421                         parameters->FindParam(&meshtype,MeshTypeEnum);
    422                         if(meshtype==Mesh2DverticalEnum){
    423                                 return CreatePVectorL2ProjectionBase();
    424                         }
    425                         else{
    426                                 return CreatePVectorL2Projection();
    427                         }
    428                         break;
    429                 #ifdef _HAVE_MASSTRANSPORT_
    430                 case MasstransportAnalysisEnum:
    431                         return CreatePVectorMasstransport();
    432                         break;
    433                 case FreeSurfaceTopAnalysisEnum:
    434                         parameters->FindParam(&meshtype,MeshTypeEnum);
    435                         switch(meshtype){
    436                                 case Mesh2DverticalEnum:
    437                                         return CreatePVectorFreeSurfaceTop1D();
    438                                 case Mesh3DEnum:
    439                                         return CreatePVectorFreeSurfaceTop();
    440                                 default:
    441                                         _error_("Mesh not supported yet");
    442                         }
    443                         break;
    444                 case FreeSurfaceBaseAnalysisEnum:
    445                         parameters->FindParam(&meshtype,MeshTypeEnum);
    446                         switch(meshtype){
    447                                 case Mesh2DverticalEnum:
    448                                         return CreatePVectorFreeSurfaceBase1D();
    449                                 case Mesh3DEnum:
    450                                         return CreatePVectorFreeSurfaceBase();
    451                                 default:
    452                                         _error_("Mesh not supported yet");
    453                         }
    454                         break;
    455                 case ExtrudeFromBaseAnalysisEnum: case ExtrudeFromTopAnalysisEnum:
    456                         return NULL;
    457                         break;
    458                 #endif
    459                 #ifdef _HAVE_HYDROLOGY_
    460                 case HydrologyShreveAnalysisEnum:
    461                         return CreatePVectorHydrologyShreve();
    462                         break;
    463                 case HydrologyDCInefficientAnalysisEnum:
    464                         return CreatePVectorHydrologyDCInefficient();
    465                         break;
    466                 case HydrologyDCEfficientAnalysisEnum:
    467                         return CreatePVectorHydrologyDCEfficient();
    468                         break;
    469           case L2ProjectionEPLAnalysisEnum:
    470                         return CreatePVectorL2ProjectionEPL();
    471                         break;
    472                 #endif
    473                 #ifdef _HAVE_BALANCED_
    474                 case BalancethicknessAnalysisEnum:
    475                         return CreatePVectorBalancethickness();
    476                         break;
    477                 case BalancevelocityAnalysisEnum:
    478                         return CreatePVectorBalancevelocity();
    479                         break;
    480                 case SmoothedSurfaceSlopeXAnalysisEnum:
    481                         return CreatePVectorSmoothedSlopeX();
    482                         break;
    483                 case SmoothedSurfaceSlopeYAnalysisEnum:
    484                         return CreatePVectorSmoothedSlopeY();
    485                         break;
    486                 #endif
    487                 #ifdef _HAVE_CONTROL_
    488                 case AdjointBalancethicknessAnalysisEnum:
    489                         return CreatePVectorAdjointBalancethickness();
    490                         break;
    491                 case AdjointHorizAnalysisEnum:
    492                         return CreatePVectorAdjointHoriz();
    493                         break;
    494                 #endif
    495                 default:
    496                         _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet");
    497         }
    498 
    499         /*make compiler happy*/
    500         return NULL;
    501 }
    502 /*}}}*/
    503 /*FUNCTION Tria::CreatePVectorL2ProjectionBase{{{*/
    504 ElementVector* Tria::CreatePVectorL2ProjectionBase(void){
    505 
    506         if(!HasEdgeOnBed()) return NULL;
    507 
    508         int index1,index2;
    509         this->EdgeOnBedIndices(&index1,&index2);
    510 
    511         Seg* seg=(Seg*)SpawnSeg(index1,index2);
    512         ElementVector* pe=seg->CreatePVectorL2Projection();
    513         delete seg->material; delete seg;
    514 
    515         /*clean up and return*/
    516         return pe;
    517 }
    518 /*}}}*/
    519 /*FUNCTION Tria::CreatePVectorL2Projection {{{*/
    520 ElementVector* Tria::CreatePVectorL2Projection(void){
    521 
    522         /*Intermediaries */
    523         int        i,input_enum;
    524         IssmDouble Jdet,value;
    525         IssmDouble xyz_list[NUMVERTICES][3];
    526         IssmDouble slopes[2];
    527         Input*     input  = NULL;
    528         Input*     input2 = NULL;
    529 
    530         /*Fetch number of nodes and dof for this finite element*/
    531         int numnodes = this->NumberofNodes();
    532 
    533         /*Initialize Element vector*/
    534         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    535         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    536 
    537         /*Retrieve all inputs and parameters*/
    538         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    539         this->parameters->FindParam(&input_enum,InputToL2ProjectEnum);
    540         switch(input_enum){
    541                 case SurfaceSlopeXEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break;
    542                 case SurfaceSlopeYEnum: input2 = inputs->GetInput(SurfaceEnum); _assert_(input2); break;
    543                 case BedSlopeXEnum:     input2 = inputs->GetInput(BedEnum);     _assert_(input2); break;
    544                 case BedSlopeYEnum:     input2 = inputs->GetInput(BedEnum);     _assert_(input2); break;
    545                 default: input = inputs->GetInput(input_enum);
    546         }
    547 
    548         /* Start  looping on the number of gaussian points: */
    549         GaussTria* gauss=new GaussTria(2);
    550         for(int ig=gauss->begin();ig<gauss->end();ig++){
    551 
    552                 gauss->GaussPoint(ig);
    553 
    554                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    555                 GetNodalFunctions(basis,gauss);
    556 
    557                 if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss);
    558                 switch(input_enum){
    559                         case SurfaceSlopeXEnum: case BedSlopeXEnum: value = slopes[0]; break;
    560                         case SurfaceSlopeYEnum: case BedSlopeYEnum: value = slopes[1]; break;
    561                         default: input->GetInputValue(&value,gauss);
    562                 }
    563 
    564                 for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
    565         }
    566 
    567         /*Clean up and return*/
    568         xDelete<IssmDouble>(basis);
    569         delete gauss;
    570         return pe;
    571 }
    572 /*}}}*/
    573376/*FUNCTION Tria::CreateJacobianMatrix{{{*/
    574377void  Tria::CreateJacobianMatrix(Matrix<IssmDouble>* Jff){
    … …  
    38173620}
    38183621/*}}}*/
    3819 /*FUNCTION Tria::CreatePVectorStressbalanceFS {{{*/
    3820 ElementVector* Tria::CreatePVectorStressbalanceFS(void){
    3821 
    3822         ElementVector* pe1;
    3823         ElementVector* pe2;
    3824         ElementVector* pe3;
    3825         ElementVector* pe;
    3826 
    3827         /*compute all stiffness matrices for this element*/
    3828         pe1=CreatePVectorStressbalanceFSViscous();
    3829         pe2=CreatePVectorStressbalanceFSShelf();
    3830         pe3=CreatePVectorStressbalanceFSFront();
    3831         pe =new ElementVector(pe1,pe2,pe3);
    3832 
    3833         /*clean-up and return*/
    3834         delete pe1;
    3835         delete pe2;
    3836         delete pe3;
    3837         return pe;
    3838 }
    3839 /*}}}*/
    3840 /*FUNCTION Tria::CreatePVectorStressbalanceFSFront{{{*/
    3841 ElementVector* Tria::CreatePVectorStressbalanceFSFront(void){
    3842 
    3843         /*If no front, return NULL*/
    3844         if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
    3845 
    3846         /*Intermediaries*/
    3847         int         i;
    3848         IssmDouble  rho_ice,rho_water,gravity,y_g;
    3849         IssmDouble  Jdet,water_pressure,air_pressure,pressure;
    3850         IssmDouble  xyz_list[NUMVERTICES][3];
    3851         IssmDouble* xyz_list_front = NULL;
    3852         IssmDouble  area_coordinates[2][3];
    3853         IssmDouble  normal[2];
    3854         GaussTria*  gauss = NULL;
    3855 
    3856         /*Fetch number of nodes and dof for this finite element*/
    3857         int vnumnodes = this->NumberofNodesVelocity();
    3858         int pnumnodes = this->NumberofNodesPressure();
    3859 
    3860         /*Prepare coordinate system list*/
    3861         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    3862         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
    3863         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    3864 
    3865         /*Initialize Element matrix and vectors*/
    3866         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3867         ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    3868         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    3869 
    3870         /*Retrieve all inputs and parameters*/
    3871         rho_water = matpar->GetRhoWater();
    3872         rho_ice = matpar->GetRhoIce();
    3873         gravity   = matpar->GetG();
    3874         ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    3875         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
    3876         NormalSection(&normal[0],xyz_list_front);
    3877 
    3878         /*Start looping on Gaussian points*/
    3879         IssmDouble ymax=max(xyz_list_front[0*3+1],xyz_list_front[1*3+1]);
    3880         IssmDouble ymin=min(xyz_list_front[0*3+1],xyz_list_front[1*3+1]);
    3881         if(ymax>0. && ymin<0.) gauss=new GaussTria(area_coordinates,30); //refined in vertical because of the sea level discontinuity
    3882         else                   gauss=new GaussTria(area_coordinates,3);
    3883 
    3884         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3885 
    3886                 gauss->GaussPoint(ig);
    3887                 y_g=GetYcoord(gauss);
    3888                 GetNodalFunctionsVelocity(vbasis,gauss);
    3889                 GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss);
    3890 
    3891                 water_pressure = rho_water*gravity*min(0.,y_g); //0 if the gaussian point is above water level
    3892                 air_pressure   = 0.;
    3893                 pressure       = water_pressure + air_pressure;
    3894 
    3895                 for(i=0;i<vnumnodes;i++){
    3896                         pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
    3897                         pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
    3898                 }
    3899         }
    3900 
    3901         /*Transform coordinate system*/
    3902         ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    3903 
    3904         /*Clean up and return*/
    3905         xDelete<int>(cs_list);
    3906         xDelete<IssmDouble>(xyz_list_front);
    3907         xDelete<IssmDouble>(vbasis);
    3908         delete gauss;
    3909         return pe;
    3910 
    3911 }
    3912 /*}}}*/
    3913 /*FUNCTION Tria::CreatePVectorStressbalanceFSViscous {{{*/
    3914 ElementVector* Tria::CreatePVectorStressbalanceFSViscous(void){
    3915 
    3916         /*Intermediaries*/
    3917         int        i;
    3918         int        approximation;
    3919         IssmDouble Jdet,gravity,rho_ice;
    3920         IssmDouble forcex,forcey;
    3921         IssmDouble xyz_list[NUMVERTICES][3];
    3922         GaussTria *gauss=NULL;
    3923 
    3924         /*Fetch number of nodes and dof for this finite element*/
    3925         int vnumnodes = this->NumberofNodesVelocity();
    3926         int pnumnodes = this->NumberofNodesPressure();
    3927 
    3928         /*Prepare coordinate system list*/
    3929         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    3930         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
    3931         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    3932 
    3933         /*Initialize Element matrix and vectors*/
    3934         ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    3935         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    3936 
    3937         /*Retrieve all inputs and parameters*/
    3938         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    3939         Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
    3940         Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
    3941         rho_ice = matpar->GetRhoIce();
    3942         gravity = matpar->GetG();
    3943 
    3944         /* Start  looping on the number of gaussian points: */
    3945         gauss=new GaussTria(5);
    3946         for(int ig=gauss->begin();ig<gauss->end();ig++){
    3947 
    3948                 gauss->GaussPoint(ig);
    3949 
    3950                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    3951                 GetNodalFunctionsVelocity(vbasis, gauss);
    3952 
    3953                 loadingforcex_input->GetInputValue(&forcex,gauss);
    3954                 loadingforcey_input->GetInputValue(&forcey,gauss);
    3955 
    3956                 for(i=0;i<vnumnodes;i++){
    3957                         pe->values[i*NDOF2+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
    3958                         pe->values[i*NDOF2+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
    3959                         pe->values[i*NDOF2+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
    3960                 }
    3961         }
    3962 
    3963         /*Transform coordinate system*/
    3964         ::TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
    3965 
    3966         /*Clean up and return*/
    3967         delete gauss;
    3968         xDelete<int>(cs_list);
    3969         xDelete<IssmDouble>(vbasis);
    3970         return pe;
    3971 }
    3972 /*}}}*/
    3973 /*FUNCTION Tria::CreatePVectorStressbalanceFSShelf{{{*/
    3974 ElementVector* Tria::CreatePVectorStressbalanceFSShelf(void){
    3975 
    3976         /*Intermediaries*/
    3977         int         i,j;
    3978         int         indices[2];
    3979         IssmDouble  gravity,rho_water,bed,water_pressure;
    3980         IssmDouble  normal_vel,vx,vy,dt;
    3981         IssmDouble      xyz_list_seg[NUMVERTICES1D][3];
    3982         IssmDouble  xyz_list[NUMVERTICES][3];
    3983         IssmDouble      normal[2];
    3984         IssmDouble  Jdet;
    3985 
    3986         /*Initialize Element vector and return if necessary*/
    3987         if(!HasEdgeOnBed() || !IsFloating()) return NULL;
    3988 
    3989         /*Fetch number of nodes and dof for this finite element*/
    3990         int vnumnodes = this->NumberofNodesVelocity();
    3991         int pnumnodes = this->NumberofNodesPressure();
    3992 
    3993         /*Prepare coordinate system list*/
    3994         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    3995         for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
    3996         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    3997 
    3998         /*Initialize Element matrix and vectors*/
    3999         ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
    4000         IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
    4001 
    4002         /*Retrieve all inputs and parameters*/
    4003         rho_water=matpar->GetRhoWater();
    4004         gravity=matpar->GetG();
    4005         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4006         Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input);
    4007 
    4008         /*Get vertex indices that lie on bed*/
    4009         this->EdgeOnBedIndices(&indices[0],&indices[1]);
    4010 
    4011         for(i=0;i<NUMVERTICES1D;i++) for(j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
    4012 
    4013         /* Start looping on the number of gauss 1d (nodes on the bedrock) */
    4014         GaussTria* gauss=new GaussTria(indices[0],indices[1],5);
    4015         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4016 
    4017                 gauss->GaussPoint(ig);
    4018 
    4019                 GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
    4020                 GetNodalFunctionsVelocity(vbasis, gauss);
    4021 
    4022                 NormalSection(&normal[0],&xyz_list_seg[0][0]);
    4023                 _assert_(normal[1]<0.);
    4024                 bed_input->GetInputValue(&bed, gauss);
    4025                 water_pressure=gravity*rho_water*bed;
    4026 
    4027                 for(i=0;i<vnumnodes;i++){
    4028                         pe->values[2*i+0]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[0];
    4029                         pe->values[2*i+1]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[1];
    4030                 }
    4031         }
    4032 
    4033         /*Transform coordinate system*/
    4034         ::TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
    4035 
    4036         /*Clean up and return*/
    4037         delete gauss;
    4038         xDelete<int>(cs_list);
    4039         xDelete<IssmDouble>(vbasis);
    4040         return pe;
    4041 }
    4042 /*}}}*/
    4043 /*FUNCTION Tria::CreatePVectorStressbalanceSSA {{{*/
    4044 ElementVector* Tria::CreatePVectorStressbalanceSSA(){
    4045 
    4046         /*compute all load vectors for this element*/
    4047         ElementVector* pe1=CreatePVectorStressbalanceSSADrivingStress();
    4048         ElementVector* pe2=CreatePVectorStressbalanceSSAFront();
    4049         ElementVector* pe =new ElementVector(pe1,pe2);
    4050 
    4051         /*clean-up and return*/
    4052         delete pe1;
    4053         delete pe2;
    4054         return pe;
    4055 }
    4056 /*}}}*/
    4057 /*FUNCTION Tria::CreatePVectorStressbalanceSSADrivingStress {{{*/
    4058 ElementVector* Tria::CreatePVectorStressbalanceSSADrivingStress(){
    4059 
    4060         /*Intermediaries */
    4061         int        i;
    4062         IssmDouble driving_stress_baseline,thickness;
    4063         IssmDouble Jdet;
    4064         IssmDouble xyz_list[NUMVERTICES][3];
    4065         IssmDouble slope[2];
    4066         IssmDouble icefrontlevel[3];
    4067 
    4068         /*Fetch number of nodes and dof for this finite element*/
    4069         int numnodes = this->NumberofNodes();
    4070 
    4071         /*Initialize Element vector and vectors*/
    4072         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,SSAApproximationEnum);
    4073         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    4074 
    4075         /*Retrieve all inputs and parameters*/
    4076         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4077         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    4078         Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
    4079         Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);
    4080         GetInputListOnVertices(&icefrontlevel[0],BedEnum);
    4081 
    4082         /* Start  looping on the number of gaussian points: */
    4083         GaussTria*     gauss  = new GaussTria(2);
    4084         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4085 
    4086                 gauss->GaussPoint(ig);
    4087 
    4088                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    4089                 GetNodalFunctions(basis, gauss);
    4090 
    4091                 thickness_input->GetInputValue(&thickness,gauss);
    4092                 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    4093                 driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;
    4094 
    4095                 /*Build load vector: */
    4096                 for(i=0;i<numnodes;i++){
    4097                         pe->values[i*NDOF2+0]+=-driving_stress_baseline*slope[0]*Jdet*gauss->weight*basis[i];
    4098                         pe->values[i*NDOF2+1]+=-driving_stress_baseline*slope[1]*Jdet*gauss->weight*basis[i];
    4099                 }
    4100         }
    4101 
    4102         /*Transform coordinate system*/
    4103         ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
    4104 
    4105         /*Clean up and return*/
    4106         xDelete<IssmDouble>(basis);
    4107         delete gauss;
    4108         return pe;
    4109 }
    4110 /*}}}*/
    4111 /*FUNCTION Tria::CreatePVectorStressbalanceSSAFront {{{*/
    4112 ElementVector* Tria::CreatePVectorStressbalanceSSAFront(){
    4113 
    4114         /*If no front, return NULL*/
    4115         if(!IsZeroLevelset(MaskIceLevelsetEnum)) return NULL;
    4116 
    4117         /*Intermediaries*/
    4118         IssmDouble  rho_ice,rho_water,gravity;
    4119         IssmDouble  Jdet,thickness,bed,water_pressure,ice_pressure,air_pressure;
    4120         IssmDouble  surface_under_water,base_under_water,pressure;
    4121         IssmDouble  xyz_list[NUMVERTICES][3];
    4122         IssmDouble  *xyz_list_front = NULL;
    4123         IssmDouble  area_coordinates[2][3];
    4124         IssmDouble  normal[2];
    4125 
    4126         /*Fetch number of nodes for this finite element*/
    4127         int numnodes = this->NumberofNodes();
    4128 
    4129         /*Initialize Element vector and other vectors*/
    4130         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    4131         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters,SSAApproximationEnum);
    4132         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    4133 
    4134         /*Retrieve all inputs and parameters*/
    4135         Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
    4136         Input* bed_input      =inputs->GetInput(BedEnum);       _assert_(bed_input);
    4137         rho_water = matpar->GetRhoWater();
    4138         rho_ice   = matpar->GetRhoIce();
    4139         gravity   = matpar->GetG();
    4140         ZeroLevelsetCoordinates(&xyz_list_front,&xyz_list[0][0],MaskIceLevelsetEnum);
    4141         GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,&xyz_list[0][0],2);
    4142         NormalSection(&normal[0],xyz_list_front);
    4143 
    4144         /*Start looping on Gaussian points*/
    4145         GaussTria* gauss=new GaussTria(area_coordinates,3);
    4146         for(int ig=gauss->begin();ig<gauss->end();ig++){
    4147 
    4148                 gauss->GaussPoint(ig);
    4149                 thickness_input->GetInputValue(&thickness,gauss);
    4150                 bed_input->GetInputValue(&bed,gauss);
    4151                 GetSegmentJacobianDeterminant(&Jdet,xyz_list_front,gauss);
    4152                 GetNodalFunctions(basis,gauss);
    4153 
    4154                 surface_under_water = min(0.,thickness+bed); // 0 if the top of the glacier is above water level
    4155                 base_under_water    = min(0.,bed);           // 0 if the bottom of the glacier is above water level
    4156                 water_pressure = 1.0/2.0*gravity*rho_water*(pow(surface_under_water,2) - pow(base_under_water,2));
    4157                 ice_pressure   = 1.0/2.0*gravity*rho_ice*pow(thickness,2);
    4158                 air_pressure   = 0;
    4159                 pressure = ice_pressure + water_pressure + air_pressure;
    4160 
    4161                 for (int i=0;i<numnodes;i++){
    4162                         pe->values[2*i+0]+= pressure*Jdet*gauss->weight*normal[0]*basis[i];
    4163                         pe->values[2*i+1]+= pressure*Jdet*gauss->weight*normal[1]*basis[i];
    4164                 }
    4165         }
    4166 
    4167         /*Transform coordinate system*/
    4168         ::TransformLoadVectorCoord(pe,nodes,numnodes,XYEnum);
    4169 
    4170         /*Clean up and return*/
    4171         xDelete<IssmDouble>(basis);
    4172         xDelete<IssmDouble>(xyz_list_front);
    4173         delete gauss;
    4174         return pe;
    4175 }
    4176 /*}}}*/
    4177 /*FUNCTION Tria::CreatePVectorStressbalanceSIA{{{*/
    4178 ElementVector* Tria::CreatePVectorStressbalanceSIA(void){
    4179 
    4180         /*Intermediaries */
    4181         IssmDouble ub,vb;
    4182         IssmDouble rho_ice,gravity,n,B,drag;
    4183         IssmDouble slope2,thickness,connectivity;
    4184         IssmDouble slope[2];
    4185 
    4186         /*Fetch number of nodes and dof for this finite element*/
    4187         int numnodes = this->NumberofNodes(); _assert_(numnodes==3);
    4188 
    4189         /*Initialize Element vector*/
    4190         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters);
    4191 
    4192         /*Retrieve all inputs and parameters*/
    4193         rho_ice=matpar->GetRhoIce();
    4194         gravity=matpar->GetG();
    4195         n=material->GetN();
    4196         B=material->GetBbar();
    4197         Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
    4198         Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
    4199         Input* thickness_input=inputs->GetInput(ThicknessEnum);  _assert_(thickness_input);
    4200         Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);  _assert_(drag_input);
    4201 
    4202         /*Spawn 3 sing elements: */
    4203         GaussTria* gauss=new GaussTria();
    4204         for(int i=0;i<numnodes;i++){
    4205 
    4206                 gauss->GaussVertex(i);
    4207 
    4208                 connectivity=(IssmDouble)vertices[i]->Connectivity();
    4209 
    4210                 thickness_input->GetInputValue(&thickness,gauss);
    4211                 drag_input->GetInputValue(&drag,gauss);
    4212                 slopex_input->GetInputValue(&slope[0],gauss);
    4213                 slopey_input->GetInputValue(&slope[1],gauss);
    4214                 slope2=slope[0]*slope[0]+slope[1]*slope[1];
    4215 
    4216                 /*Payne 1995 (m = 1, B = 5e-3/yts = 1.58e-10  m/s/Pa*/
    4217                 ub=-1.58*1.e-10*rho_ice*gravity*thickness*slope[0];
    4218                 vb=-1.58*1.e-10*rho_ice*gravity*thickness*slope[1];
    4219                 ///*Ritz et al. 1996*/
    4220                 //ub=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[0]/sqrt(slope2);
    4221                 //vb=drag*(rho_ice*gravity*thickness)*(rho_ice*gravity*thickness)*slope[1]/sqrt(slope2);
    4222                 ///*Rutt et al. 2009*/
    4223                 //ub=-drag*rho_ice*gravity*thickness*slope[0];
    4224                 //vb=-drag*rho_ice*gravity*thickness*slope[1];
    4225 
    4226                 pe->values[2*i+0]=(ub-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/connectivity;
    4227                 pe->values[2*i+1]=(vb-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/connectivity;
    4228         }
    4229 
    4230         /*Clean up and return*/
    4231         delete gauss;
    4232         return pe;
    4233 }
    4234 /*}}}*/
    42353622/*FUNCTION Tria::CreateJacobianStressbalanceSSA{{{*/
    42363623ElementMatrix* Tria::CreateJacobianStressbalanceSSA(void){
    … …  
    56074994}
    56084995/*}}}*/
    5609 /*FUNCTION Tria::CreatePVectorAdjointBalancethickness{{{*/
    5610 ElementVector* Tria::CreatePVectorAdjointBalancethickness(void){
    5611 
    5612         /*Intermediaries */
    5613         int        i,resp;
    5614         IssmDouble Jdet;
    5615         IssmDouble thickness,thicknessobs,weight;
    5616         int        num_responses;
    5617         IssmDouble xyz_list[NUMVERTICES][3];
    5618         IssmDouble dH[2];
    5619         IssmDouble vx,vy,vel;
    5620         int       *responses = NULL;
    5621 
    5622         /*Fetch number of nodes and dof for this finite element*/
    5623         int numnodes = this->NumberofNodes();
    5624 
    5625         /*Initialize Element vector and vectors*/
    5626         ElementVector* pe     = new ElementVector(nodes,numnodes,this->parameters);
    5627         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    5628         IssmDouble*    dbasis = xNew<IssmDouble>(2*numnodes);
    5629 
    5630         /*Retrieve all inputs and parameters*/
    5631         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5632         this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    5633         this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    5634         Input* thickness_input    = inputs->GetInput(ThicknessEnum);                          _assert_(thickness_input);
    5635         Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);              _assert_(thicknessobs_input);
    5636         Input* weights_input      = inputs->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input);
    5637         Input* vx_input           = inputs->GetInput(VxEnum);                                 _assert_(vx_input);
    5638         Input* vy_input           = inputs->GetInput(VyEnum);                                 _assert_(vy_input);
    5639 
    5640         /* Start  looping on the number of gaussian points: */
    5641         GaussTria* gauss=new GaussTria(2);
    5642         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5643 
    5644                 gauss->GaussPoint(ig);
    5645 
    5646                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5647                 GetNodalFunctions(basis, gauss);
    5648                 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
    5649 
    5650                 thickness_input->GetInputValue(&thickness, gauss);
    5651                 thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
    5652                 thicknessobs_input->GetInputValue(&thicknessobs, gauss);
    5653 
    5654                 /*Loop over all requested responses*/
    5655                 for(resp=0;resp<num_responses;resp++){
    5656 
    5657                         weights_input->GetInputValue(&weight,gauss,responses[resp]);
    5658 
    5659                         switch(responses[resp]){
    5660                                 case ThicknessAbsMisfitEnum:
    5661                                         for(i=0;i<numnodes;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
    5662                                         break;
    5663                                 case ThicknessAbsGradientEnum:
    5664                                         for(i=0;i<numnodes;i++) pe->values[i]+= - weight*dH[0]*dbasis[0*numnodes+i]*Jdet*gauss->weight;
    5665                                         for(i=0;i<numnodes;i++) pe->values[i]+= - weight*dH[1]*dbasis[1*numnodes+i]*Jdet*gauss->weight;
    5666                                         break;
    5667                                 case ThicknessAlongGradientEnum:
    5668                                         vx_input->GetInputValue(&vx,gauss);
    5669                                         vy_input->GetInputValue(&vy,gauss);
    5670                                         vel = sqrt(vx*vx+vy*vy);
    5671                                         vx  = vx/(vel+1.e-9);
    5672                                         vy  = vy/(vel+1.e-9);
    5673                                         for(i=0;i<numnodes;i++) pe->values[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numnodes+i]*vx+dbasis[1*numnodes+i]*vy)*Jdet*gauss->weight;
    5674                                         break;
    5675                                 case ThicknessAcrossGradientEnum:
    5676                                         vx_input->GetInputValue(&vx,gauss);
    5677                                         vy_input->GetInputValue(&vy,gauss);
    5678                                         vel = sqrt(vx*vx+vy*vy);
    5679                                         vx  = vx/(vel+1.e-9);
    5680                                         vy  = vy/(vel+1.e-9);
    5681                                         for(i=0;i<numnodes;i++) pe->values[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numnodes+i]*(-vy)+dbasis[1*numnodes+i]*vx)*Jdet*gauss->weight;
    5682                                         break;
    5683                                 default:
    5684                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    5685                         }
    5686                 }
    5687         }
    5688 
    5689         /*Clean up and return*/
    5690         xDelete<IssmDouble>(basis);
    5691         xDelete<IssmDouble>(dbasis);
    5692         xDelete<int>(responses);
    5693         delete gauss;
    5694         return pe;
    5695 }
    5696 /*}}}*/
    5697 /*FUNCTION Tria::CreatePVectorAdjointHoriz{{{*/
    5698 ElementVector* Tria::CreatePVectorAdjointHoriz(void){
    5699 
    5700         /*Intermediaries */
    5701         int        i,resp;
    5702         int       *responses=NULL;
    5703         int        num_responses;
    5704         IssmDouble Jdet;
    5705         IssmDouble obs_velocity_mag,velocity_mag;
    5706         IssmDouble dux,duy;
    5707         IssmDouble epsvel=2.220446049250313e-16;
    5708         IssmDouble meanvel=3.170979198376458e-05; /*1000 m/yr*/
    5709         IssmDouble scalex=0.,scaley=0.,scale=0.,S=0.;
    5710         IssmDouble vx,vy,vxobs,vyobs,weight;
    5711         IssmDouble xyz_list[NUMVERTICES][3];
    5712 
    5713         /*Fetch number of nodes and dof for this finite element*/
    5714         int numnodes = this->NumberofNodes();
    5715 
    5716         /*Initialize Element vector*/
    5717         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    5718         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    5719 
    5720         /*Retrieve all inputs and parameters*/
    5721         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    5722         this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
    5723         this->parameters->FindParam(&responses,NULL,InversionCostFunctionsEnum);
    5724         Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
    5725         Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
    5726         Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
    5727         Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
    5728         Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
    5729 
    5730         /*Get Surface if required by one response*/
    5731         for(resp=0;resp<num_responses;resp++){
    5732                 if(responses[resp]==SurfaceAverageVelMisfitEnum){
    5733                         inputs->GetInputValue(&S,SurfaceAreaEnum); break;
    5734                 }
    5735         }
    5736 
    5737         /* Start  looping on the number of gaussian points: */
    5738         GaussTria* gauss=new GaussTria(4);
    5739         for(int ig=gauss->begin();ig<gauss->end();ig++){
    5740 
    5741                 gauss->GaussPoint(ig);
    5742 
    5743                 /* Get Jacobian determinant: */
    5744                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    5745 
    5746                 /*Get all parameters at gaussian point*/
    5747                 vx_input->GetInputValue(&vx,gauss);
    5748                 vy_input->GetInputValue(&vy,gauss);
    5749                 vxobs_input->GetInputValue(&vxobs,gauss);
    5750                 vyobs_input->GetInputValue(&vyobs,gauss);
    5751                 GetNodalFunctions(basis, gauss);
    5752 
    5753                 /*Loop over all requested responses*/
    5754                 for(resp=0;resp<num_responses;resp++){
    5755 
    5756                         weights_input->GetInputValue(&weight,gauss,responses[resp]);
    5757 
    5758                         switch(responses[resp]){
    5759                                 case SurfaceAbsVelMisfitEnum:
    5760                                         /*
    5761                                          *      1  [           2              2 ]
    5762                                          * J = --- | (u - u   )  +  (v - v   )  |
    5763                                          *      2  [       obs            obs   ]
    5764                                          *
    5765                                          *        dJ
    5766                                          * DU = - -- = (u   - u )
    5767                                          *        du     obs
    5768                                          */
    5769                                         for(i=0;i<numnodes;i++){
    5770                                                 dux=vxobs-vx;
    5771                                                 duy=vyobs-vy;
    5772                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
    5773                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    5774                                         }
    5775                                         break;
    5776                                 case SurfaceRelVelMisfitEnum:
    5777                                         /*
    5778                                          *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
    5779                                          * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
    5780                                          *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
    5781                                          *              obs                        obs                     
    5782                                          *
    5783                                          *        dJ     \bar{v}^2
    5784                                          * DU = - -- = ------------- (u   - u )
    5785                                          *        du   (u   + eps)^2    obs
    5786                                          *               obs
    5787                                          */
    5788                                         for(i=0;i<numnodes;i++){
    5789                                                 scalex=pow(meanvel/(vxobs+epsvel),2); if(vxobs==0)scalex=0;
    5790                                                 scaley=pow(meanvel/(vyobs+epsvel),2); if(vyobs==0)scaley=0;
    5791                                                 dux=scalex*(vxobs-vx);
    5792                                                 duy=scaley*(vyobs-vy);
    5793                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
    5794                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    5795                                         }
    5796                                         break;
    5797                                 case SurfaceLogVelMisfitEnum:
    5798                                         /*
    5799                                          *                 [        vel + eps     ] 2
    5800                                          * J = 4 \bar{v}^2 | log ( -----------  ) | 
    5801                                          *                 [       vel   + eps    ]
    5802                                          *                            obs
    5803                                          *
    5804                                          *        dJ                 2 * log(...)
    5805                                          * DU = - -- = - 4 \bar{v}^2 -------------  u
    5806                                          *        du                 vel^2 + eps
    5807                                          *           
    5808                                          */
    5809                                         for(i=0;i<numnodes;i++){
    5810                                                 velocity_mag    =sqrt(pow(vx,   2)+pow(vy,   2))+epsvel;
    5811                                                 obs_velocity_mag=sqrt(pow(vxobs,2)+pow(vyobs,2))+epsvel;
    5812                                                 scale=-8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
    5813                                                 dux=scale*vx;
    5814                                                 duy=scale*vy;
    5815                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
    5816                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    5817                                         }
    5818                                         break;
    5819                                 case SurfaceAverageVelMisfitEnum:
    5820                                         /*
    5821                                          *      1                    2              2
    5822                                          * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
    5823                                          *      S                obs            obs
    5824                                          *
    5825                                          *        dJ      1       1
    5826                                          * DU = - -- = - --- ----------- * 2 (u - u   )
    5827                                          *        du      S  2 sqrt(...)           obs
    5828                                          */
    5829                                         for(i=0;i<numnodes;i++){
    5830                                                 scale=1./(S*2*sqrt(pow(vx-vxobs,2)+pow(vy-vyobs,2))+epsvel);
    5831                                                 dux=scale*(vxobs-vx);
    5832                                                 duy=scale*(vyobs-vy);
    5833                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
    5834                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    5835                                         }
    5836                                         break;
    5837                                 case SurfaceLogVxVyMisfitEnum:
    5838                                         /*
    5839                                          *      1            [        |u| + eps     2          |v| + eps     2  ]
    5840                                          * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   | 
    5841                                          *      2            [       |u    |+ eps              |v    |+ eps     ]
    5842                                          *                              obs                       obs
    5843                                          *        dJ                              1      u                             1
    5844                                          * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
    5845                                          *        du                         |u| + eps  |u|                           u + eps
    5846                                          */
    5847                                         for(i=0;i<numnodes;i++){
    5848                                                 dux = - meanvel*meanvel * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
    5849                                                 duy = - meanvel*meanvel * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
    5850                                                 pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i];
    5851                                                 pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i];
    5852                                         }
    5853                                         break;
    5854                                 case DragCoefficientAbsGradientEnum:
    5855                                         /*Nothing in P vector*/
    5856                                         break;
    5857                                 case ThicknessAbsGradientEnum:
    5858                                         /*Nothing in P vector*/
    5859                                         break;
    5860                                 case ThicknessAlongGradientEnum:
    5861                                         /*Nothing in P vector*/
    5862                                         break;
    5863                                 case ThicknessAcrossGradientEnum:
    5864                                         /*Nothing in P vector*/
    5865                                         break;
    5866                                 case RheologyBbarAbsGradientEnum:
    5867                                         /*Nothing in P vector*/
    5868                                         break;
    5869                                 default:
    5870                                         _error_("response " << EnumToStringx(responses[resp]) << " not supported yet");
    5871                         }
    5872                 }
    5873         }
    5874 
    5875         /*Clean up and return*/
    5876         xDelete<IssmDouble>(basis);
    5877         xDelete<int>(responses);
    5878         delete gauss;
    5879         return pe;
    5880 }
    5881 /*}}}*/
    58824996/*FUNCTION Tria::DragCoefficientAbsGradient{{{*/
    58834997IssmDouble Tria::DragCoefficientAbsGradient(void){
    … …  
    64495563        xDelete<IssmDouble>(basis);
    64505564        return Ke;
    6451 }
    6452 /*}}}*/
    6453 /*FUNCTION Tria::CreatePVectorHydrologyShreve {{{*/
    6454 ElementVector* Tria::CreatePVectorHydrologyShreve(void){
    6455 
    6456         /*Intermediaries */
    6457         int        i;
    6458         IssmDouble Jdet,dt;
    6459         IssmDouble basal_melting_g;
    6460         IssmDouble old_watercolumn_g;
    6461         IssmDouble xyz_list[NUMVERTICES][3];
    6462 
    6463         /*Skip if water or ice shelf element*/
    6464         if(NoIceInElement() || IsFloating()) return NULL;
    6465 
    6466         /*Fetch number of nodes and dof for this finite element*/
    6467         int numnodes = this->NumberofNodes();
    6468 
    6469         /*Initialize Element vector*/
    6470         ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters);
    6471         IssmDouble*    basis = xNewZeroInit<IssmDouble>(numnodes);
    6472 
    6473         /*Retrieve all inputs and parameters*/
    6474         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6475         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6476         Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
    6477         Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);         _assert_(old_watercolumn_input);
    6478         /*Initialize basal_melting_correction_g to 0, do not forget!:*/
    6479         /* Start  looping on the number of gaussian points: */
    6480         GaussTria* gauss=new GaussTria(2);
    6481         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6482 
    6483                 gauss->GaussPoint(ig);
    6484 
    6485                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6486                 GetNodalFunctions(basis, gauss);
    6487 
    6488                 basal_melting_input->GetInputValue(&basal_melting_g,gauss);
    6489                 old_watercolumn_input->GetInputValue(&old_watercolumn_g,gauss);
    6490 
    6491                 if(reCast<int,IssmDouble>(dt)){
    6492                         for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
    6493                 }
    6494                 else{
    6495                         for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*basal_melting_g*basis[i];
    6496                 }
    6497         }
    6498 
    6499         /*Clean up and return*/
    6500         xDelete<IssmDouble>(basis);
    6501         delete gauss;
    6502         return pe;
    6503 }
    6504 /*}}}*/
    6505 /*FUNCTION Tria::CreatePVectorHydrologyDCInefficient {{{*/
    6506 ElementVector* Tria::CreatePVectorHydrologyDCInefficient(void){
    6507 
    6508         /*Intermediaries */
    6509         IssmDouble Jdet;
    6510         IssmDouble xyz_list[NUMVERTICES][3];
    6511         IssmDouble dt,scalar,water_head;
    6512         IssmDouble water_load,transfer;
    6513         IssmDouble sediment_storing;
    6514 
    6515         /*Fetch number of nodes and dof for this finite element*/
    6516         int numnodes = this->NumberofNodes();
    6517 
    6518         /*Initialize Element vector and other vectors*/
    6519         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    6520         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    6521 
    6522         /*Retrieve all inputs and parameters*/
    6523         sediment_storing = matpar->GetSedimentStoring();
    6524         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6525         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6526         Input* water_input=inputs->GetInput(BasalforcingsMeltingRateEnum);  _assert_(water_input);
    6527         Input* transfer_input=inputs->GetInput(WaterTransferEnum);  _assert_(transfer_input);
    6528         Input* old_wh_input=NULL;
    6529 
    6530         if(reCast<bool,IssmDouble>(dt)){
    6531                 old_wh_input=inputs->GetInput(SedimentHeadOldEnum); _assert_(old_wh_input);
    6532         }
    6533 
    6534         /* Start  looping on the number of gaussian points: */
    6535         GaussTria* gauss=new GaussTria(2);
    6536         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6537 
    6538                 gauss->GaussPoint(ig);
    6539                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6540                 GetNodalFunctions(basis, gauss);
    6541 
    6542                 /*Loading term*/
    6543                 water_input->GetInputValue(&water_load,gauss);
    6544                 transfer_input->GetInputValue(&transfer,gauss);
    6545                 scalar = Jdet*gauss->weight*(water_load+transfer);
    6546                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
    6547                 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    6548 
    6549                 /*Transient term*/
    6550                 if(reCast<bool,IssmDouble>(dt)){
    6551                         old_wh_input->GetInputValue(&water_head,gauss);
    6552                         scalar = Jdet*gauss->weight*water_head*sediment_storing;
    6553                         for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    6554                 }
    6555         }
    6556 
    6557         /*Clean up and return*/
    6558         xDelete<IssmDouble>(basis);
    6559         delete gauss;
    6560         return pe;
    6561 }
    6562 /*}}}*/
    6563 /*FUNCTION Tria::CreatePVectorHydrologyDCEfficient {{{*/
    6564 ElementVector* Tria::CreatePVectorHydrologyDCEfficient(void){
    6565 
    6566         /*Intermediaries */
    6567         IssmDouble connectivity;
    6568         IssmDouble Jdet;
    6569         IssmDouble xyz_list[NUMVERTICES][3];
    6570         IssmDouble dt,scalar,water_head;
    6571         IssmDouble transfer,residual;
    6572         IssmDouble epl_thickness;
    6573         IssmDouble epl_specificstoring;
    6574         GaussTria* gauss = NULL;
    6575 
    6576         /*Check that all nodes are active, else return empty matrix*/
    6577         if(!this->AllActive()){
    6578                 return NULL;
    6579         }
    6580 
    6581         /*Fetch number of nodes and dof for this finite element*/
    6582         int numnodes = this->NumberofNodes();
    6583 
    6584         /*Initialize Element vector and other vectors*/
    6585         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    6586         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    6587 
    6588         /*Retrieve all inputs and parameters*/
    6589         epl_specificstoring = matpar->GetEplSpecificStoring();
    6590         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6591         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    6592         Input* residual_input=inputs->GetInput(SedimentHeadResidualEnum);     _assert_(residual_input);
    6593         Input* transfer_input=inputs->GetInput(WaterTransferEnum);            _assert_(transfer_input);
    6594         Input* thickness_input=inputs->GetInput(HydrologydcEplThicknessEnum); _assert_(thickness_input);
    6595         Input* old_wh_input=NULL;
    6596 
    6597         if(reCast<bool,IssmDouble>(dt)){
    6598                 old_wh_input=inputs->GetInput(EplHeadOldEnum); _assert_(old_wh_input);
    6599         }
    6600 
    6601         /* Start  looping on the number of gaussian points: */
    6602         gauss=new GaussTria(2);
    6603         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6604 
    6605                 gauss->GaussPoint(ig);
    6606                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6607                 GetNodalFunctions(basis,gauss);
    6608 
    6609                 /*Loading term*/
    6610                 transfer_input->GetInputValue(&transfer,gauss);
    6611                 scalar = Jdet*gauss->weight*(-transfer);
    6612                 if(reCast<bool,IssmDouble>(dt)) scalar = scalar*dt;
    6613                 for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    6614 
    6615                 /*Transient term*/
    6616                 if(reCast<bool,IssmDouble>(dt)){
    6617                         thickness_input->GetInputValue(&epl_thickness,gauss);
    6618                         old_wh_input->GetInputValue(&water_head,gauss);
    6619                         scalar = Jdet*gauss->weight*water_head*epl_specificstoring*epl_thickness;
    6620                         for(int i=0;i<numnodes;i++) pe->values[i]+=scalar*basis[i];
    6621                 }
    6622         }
    6623         delete gauss;
    6624 
    6625         /*      Add residual if necessary*/
    6626         gauss=new GaussTria();
    6627         for(int iv=0;iv<NUMVERTICES;iv++){
    6628                 gauss->GaussVertex(iv);
    6629 
    6630                 connectivity = IssmDouble(vertices[iv]->Connectivity());
    6631                 residual_input->GetInputValue(&residual,gauss);
    6632                 pe->values[iv]+=residual/connectivity;
    6633         }
    6634 
    6635         /*Clean up and return*/
    6636         xDelete<IssmDouble>(basis);
    6637         delete gauss;
    6638         return pe;
    6639 }
    6640 /*}}}*/
    6641 /*FUNCTION Tria::CreatePVectorL2ProjectionEPL {{{*/
    6642 ElementVector* Tria::CreatePVectorL2ProjectionEPL(void){
    6643 
    6644         /*Intermediaries */
    6645         int        i,input_enum;
    6646         IssmDouble Jdet,value;
    6647         IssmDouble xyz_list[NUMVERTICES][3];
    6648         IssmDouble slopes[2];
    6649         Input*     input  = NULL;
    6650         Input*     input2 = NULL;
    6651 
    6652         /*Fetch number of nodes and dof for this finite element*/
    6653         int numnodes = this->NumberofNodes();
    6654 
    6655         /*Initialize Element vector*/
    6656         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    6657         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    6658 
    6659         /*Retrieve all inputs and parameters*/
    6660         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    6661         this->parameters->FindParam(&input_enum,InputToL2ProjectEnum);
    6662         switch(input_enum){
    6663                 case EplHeadSlopeXEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break;
    6664                 case EplHeadSlopeYEnum: input2 = inputs->GetInput(EplHeadEnum); _assert_(input2); break;
    6665         default: input = inputs->GetInput(input_enum);
    6666         }
    6667 
    6668         /* Start  looping on the number of gaussian points: */
    6669         GaussTria* gauss=new GaussTria(2);
    6670         for(int ig=gauss->begin();ig<gauss->end();ig++){
    6671 
    6672                 gauss->GaussPoint(ig);
    6673 
    6674                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    6675                 GetNodalFunctions(basis,gauss);
    6676 
    6677                 if(input2) input2->GetInputDerivativeValue(&slopes[0],&xyz_list[0][0],gauss);
    6678                 switch(input_enum){
    6679                         case EplHeadSlopeXEnum: value = slopes[0]; break;
    6680                         case EplHeadSlopeYEnum: value = slopes[1]; break;
    6681                         default: input->GetInputValue(&value,gauss);
    6682                 }
    6683 
    6684                 for(i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*value*basis[i];
    6685         }
    6686 
    6687         /*Clean up and return*/
    6688         xDelete<IssmDouble>(basis);
    6689         delete gauss;
    6690         return pe;
    66915565}
    66925566/*}}}*/
    … …  
    76166490}
    76176491/*}}}*/
    7618 /*FUNCTION Tria::CreatePVectorMasstransport{{{*/
    7619 ElementVector* Tria::CreatePVectorMasstransport(void){
    7620 
    7621         switch(GetElementType()){
    7622                 case P1Enum: case P2Enum:
    7623                         return CreatePVectorMasstransport_CG();
    7624                 case P1DGEnum:
    7625                         return CreatePVectorMasstransport_DG();
    7626                 default:
    7627                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    7628         }
    7629 }
    7630 /*}}}*/
    7631 /*FUNCTION Tria::CreatePVectorMasstransport_CG {{{*/
    7632 ElementVector* Tria::CreatePVectorMasstransport_CG(void){
    7633 
    7634         /*Intermediaries */
    7635         IssmDouble Jdet,dt;
    7636         IssmDouble ms,mb,mb_correction,thickness;
    7637         IssmDouble xyz_list[NUMVERTICES][3];
    7638 
    7639         /*Fetch number of nodes and dof for this finite element*/
    7640         int numnodes = this->NumberofNodes();
    7641 
    7642         /*Initialize Element vector and other vectors*/
    7643         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    7644         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    7645 
    7646         /*Retrieve all inputs and parameters*/
    7647         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7648         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    7649         Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum);      _assert_(ms_input);
    7650         Input* mb_input     = inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(mb_input);
    7651         Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
    7652         Input* thickness_input  = inputs->GetInput(ThicknessEnum);     _assert_(thickness_input);
    7653 
    7654         /*Initialize mb_correction to 0, do not forget!:*/
    7655         /* Start  looping on the number of gaussian points: */
    7656         GaussTria* gauss=new GaussTria(2);
    7657         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7658 
    7659                 gauss->GaussPoint(ig);
    7660 
    7661                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7662                 GetNodalFunctions(basis,gauss);
    7663 
    7664                 ms_input->GetInputValue(&ms,gauss);
    7665                 mb_input->GetInputValue(&mb,gauss);
    7666                 thickness_input->GetInputValue(&thickness,gauss);
    7667                 if(mb_correction_input)
    7668                  mb_correction_input->GetInputValue(&mb_correction,gauss);
    7669                 else
    7670                  mb_correction=0.;
    7671 
    7672                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb-mb_correction))*basis[i];
    7673         }
    7674 
    7675         /*Clean up and return*/
    7676         xDelete<IssmDouble>(basis);
    7677         delete gauss;
    7678         return pe;
    7679 }
    7680 /*}}}*/
    7681 /*FUNCTION Tria::CreatePVectorMasstransport_DG {{{*/
    7682 ElementVector* Tria::CreatePVectorMasstransport_DG(void){
    7683 
    7684         /*Intermediaries */
    7685         IssmDouble Jdet,dt;
    7686         IssmDouble ms,mb,thickness;
    7687         IssmDouble xyz_list[NUMVERTICES][3];
    7688 
    7689         /*Fetch number of nodes and dof for this finite element*/
    7690         int numnodes = this->NumberofNodes();
    7691 
    7692         /*Initialize Element vector and other vectors*/
    7693         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    7694         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    7695 
    7696         /*Retrieve all inputs and parameters*/
    7697         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    7698         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7699         Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
    7700         Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
    7701         Input* thickness_input=inputs->GetInput(ThicknessEnum);           _assert_(thickness_input);
    7702 
    7703         /* Start  looping on the number of gaussian points: */
    7704         GaussTria* gauss=new GaussTria(2);
    7705         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7706 
    7707                 gauss->GaussPoint(ig);
    7708 
    7709                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7710                 GetNodalFunctions(basis,gauss);
    7711 
    7712                 ms_input->GetInputValue(&ms,gauss);
    7713                 mb_input->GetInputValue(&mb,gauss);
    7714                 thickness_input->GetInputValue(&thickness,gauss);
    7715 
    7716                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(thickness+dt*(ms-mb))*basis[i];
    7717         }
    7718 
    7719         /*Clean up and return*/
    7720         xDelete<IssmDouble>(basis);
    7721         delete gauss;
    7722         return pe;
    7723 }
    7724 /*}}}*/
    7725 /*FUNCTION Tria::CreatePVectorFreeSurfaceTop {{{*/
    7726 ElementVector* Tria::CreatePVectorFreeSurfaceTop(void){
    7727 
    7728         /*Intermediaries */
    7729         IssmDouble Jdet,dt;
    7730         IssmDouble ms,surface,vz;
    7731         IssmDouble xyz_list[NUMVERTICES][3];
    7732 
    7733         /*Fetch number of nodes and dof for this finite element*/
    7734         int numnodes = this->NumberofNodes();
    7735 
    7736         /*Initialize Element vector and other vectors*/
    7737         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    7738         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    7739 
    7740         /*Retrieve all inputs and parameters*/
    7741         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7742         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    7743         Input* vz_input     = inputs->GetInput(VzEnum);                         _assert_(vz_input);
    7744         Input* ms_input     = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
    7745         Input* surface_input= inputs->GetInput(SurfaceEnum);                    _assert_(surface_input);
    7746 
    7747         /*Initialize mb_correction to 0, do not forget!:*/
    7748         /* Start  looping on the number of gaussian points: */
    7749         GaussTria* gauss=new GaussTria(2);
    7750         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7751 
    7752                 gauss->GaussPoint(ig);
    7753 
    7754                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7755                 GetNodalFunctions(basis,gauss);
    7756 
    7757                 vz_input->GetInputValue(&vz,gauss);
    7758                 ms_input->GetInputValue(&ms,gauss);
    7759                 surface_input->GetInputValue(&surface,gauss);
    7760 
    7761                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(surface + dt*ms + dt*vz)*basis[i];
    7762         }
    7763 
    7764         /*Clean up and return*/
    7765         xDelete<IssmDouble>(basis);
    7766         delete gauss;
    7767         return pe;
    7768 }
    7769 /*}}}*/
    7770 /*FUNCTION Tria::CreatePVectorFreeSurfaceTop1D {{{*/
    7771 ElementVector* Tria::CreatePVectorFreeSurfaceTop1D(void){
    7772 
    7773         if(!HasEdgeOnSurface()) return NULL;
    7774 
    7775         int index1,index2;
    7776         this->EdgeOnSurfaceIndices(&index1,&index2);
    7777 
    7778         Seg* seg=(Seg*)SpawnSeg(index1,index2);
    7779         ElementVector* pe=seg->CreatePVectorFreeSurfaceTop();
    7780         delete seg->material; delete seg;
    7781 
    7782         /*clean up and return*/
    7783         return pe;
    7784 }
    7785 /*}}}*/
    7786 /*FUNCTION Tria::CreatePVectorFreeSurfaceBase {{{*/
    7787 ElementVector* Tria::CreatePVectorFreeSurfaceBase(void){
    7788 
    7789         /*Intermediaries */
    7790         IssmDouble Jdet,dt;
    7791         IssmDouble mb,mb_correction,bed,vz;
    7792         IssmDouble xyz_list[NUMVERTICES][3];
    7793 
    7794         /*Fetch number of nodes and dof for this finite element*/
    7795         int numnodes = this->NumberofNodes();
    7796 
    7797         /*Initialize Element vector and other vectors*/
    7798         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    7799         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    7800 
    7801         /*Retrieve all inputs and parameters*/
    7802         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    7803         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    7804         Input* vz_input  = inputs->GetInput(VzEnum);                         _assert_(vz_input);
    7805         Input* mb_input  = inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
    7806         Input* mb_correction_input = inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
    7807         Input* bed_input = inputs->GetInput(BedEnum);                        _assert_(bed_input);
    7808 
    7809         /*Initialize mb_correction to 0, do not forget!:*/
    7810         /* Start  looping on the number of gaussian points: */
    7811         GaussTria* gauss=new GaussTria(2);
    7812         for(int ig=gauss->begin();ig<gauss->end();ig++){
    7813 
    7814                 gauss->GaussPoint(ig);
    7815 
    7816                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    7817                 GetNodalFunctions(basis,gauss);
    7818 
    7819                 vz_input->GetInputValue(&vz,gauss);
    7820                 mb_input->GetInputValue(&mb,gauss);
    7821                 bed_input->GetInputValue(&bed,gauss);
    7822                 if(mb_correction_input)
    7823                  mb_correction_input->GetInputValue(&mb_correction,gauss);
    7824                 else
    7825                  mb_correction=0.;
    7826 
    7827                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(bed+dt*(mb-mb_correction) + dt*vz)*basis[i];
    7828         }
    7829 
    7830         /*Clean up and return*/
    7831         xDelete<IssmDouble>(basis);
    7832         delete gauss;
    7833         return pe;
    7834 }
    7835 /*}}}*/
    7836 /*FUNCTION Tria::CreatePVectorFreeSurfaceBase1D {{{*/
    7837 ElementVector* Tria::CreatePVectorFreeSurfaceBase1D(void){
    7838 
    7839         if(!HasEdgeOnBed()) return NULL;
    7840 
    7841         int index1,index2;
    7842         this->EdgeOnBedIndices(&index1,&index2);
    7843 
    7844         Seg* seg=(Seg*)SpawnSeg(index1,index2);
    7845         ElementVector* pe=seg->CreatePVectorFreeSurfaceBase();
    7846         delete seg->material; delete seg;
    7847 
    7848         /*clean up and return*/
    7849         return pe;
    7850 }
    7851 /*}}}*/
    78526492#endif
    78536493
    … …  
    79826622        delete gauss;
    79836623        return Ke;
    7984 }
    7985 /*}}}*/
    7986 /*FUNCTION Tria::CreatePVectorDamageEvolution{{{*/
    7987 ElementVector* Tria::CreatePVectorDamageEvolution(void){
    7988 
    7989         switch(GetElementType()){
    7990                 case P1Enum: case P2Enum:
    7991                         return CreatePVectorDamageEvolution_CG();
    7992                 case P1DGEnum:
    7993                         _error_("DG not implemented yet");
    7994                 default:
    7995                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    7996         }
    7997 }
    7998 /*}}}*/
    7999 /*FUNCTION Tria::CreatePVectorDamageEvolution_CG {{{*/
    8000 ElementVector* Tria::CreatePVectorDamageEvolution_CG(void){
    8001 
    8002         /*Intermediaries */
    8003         IssmDouble  Jdet ,dt;
    8004         IssmDouble  f,damage;
    8005         IssmDouble  xyz_list[NUMVERTICES][3];
    8006         IssmDouble  f_list[NUMVERTICES];
    8007         Input      *damage_input             = NULL;
    8008 
    8009         /*Fetch number of nodes and dof for this finite element*/
    8010         int numnodes = this->NumberofNodes();
    8011 
    8012         /*Initialize Element vector and other vectors*/
    8013         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    8014         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    8015 
    8016         /*Retrieve all inputs and parameters*/
    8017         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8018         this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    8019         damage_input  = this->material->inputs->GetInput(DamageDbarEnum);     _assert_(damage_input);
    8020        
    8021         /*retrieve damage evolution forcing function: */
    8022         this->DamageEvolutionF(&f_list[0]);
    8023 
    8024         /*Initialize forcing function f to 0, do not forget!:*/
    8025         /* Start  looping on the number of gaussian points: */
    8026         GaussTria* gauss=new GaussTria(2);
    8027         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8028 
    8029                 gauss->GaussPoint(ig);
    8030 
    8031                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8032                 GetNodalFunctions(basis,gauss);
    8033 
    8034                 TriaRef::GetInputValue(&f,f_list,gauss);
    8035                 damage_input->GetInputValue(&damage,gauss);
    8036 
    8037                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
    8038         }
    8039 
    8040         /*Clean up and return*/
    8041         xDelete<IssmDouble>(basis);
    8042         delete gauss;
    8043         return pe;
    80446624}
    80456625/*}}}*/
    … …  
    85507130        xDelete<IssmDouble>(basis);
    85517131        return Ke;
    8552 }
    8553 /*}}}*/
    8554 /*FUNCTION Tria::CreatePVectorBalancethickness{{{*/
    8555 ElementVector* Tria::CreatePVectorBalancethickness(void){
    8556 
    8557         switch(GetElementType()){
    8558                 case P1Enum:
    8559                         return CreatePVectorBalancethickness_CG();
    8560                         break;
    8561                 case P1DGEnum:
    8562                         return CreatePVectorBalancethickness_DG();
    8563                 default:
    8564                         _error_("Element type " << EnumToStringx(GetElementType()) << " not supported yet");
    8565         }
    8566 }
    8567 /*}}}*/
    8568 /*FUNCTION Tria::CreatePVectorBalancethickness_CG{{{*/
    8569 ElementVector* Tria::CreatePVectorBalancethickness_CG(void){
    8570 
    8571         /*Intermediaries */
    8572         IssmDouble xyz_list[NUMVERTICES][3];
    8573         IssmDouble dhdt_g,mb_g,ms_g,Jdet;
    8574 
    8575         /*Fetch number of nodes and dof for this finite element*/
    8576         int numnodes = this->NumberofNodes();
    8577 
    8578         /*Initialize Element vector and other vectors*/
    8579         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    8580         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    8581 
    8582         /*Retrieve all inputs and parameters*/
    8583         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8584         Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
    8585         Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
    8586         Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
    8587 
    8588         /* Start  looping on the number of gaussian points: */
    8589         GaussTria* gauss=new GaussTria(2);
    8590         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8591 
    8592                 gauss->GaussPoint(ig);
    8593 
    8594                 ms_input->GetInputValue(&ms_g,gauss);
    8595                 mb_input->GetInputValue(&mb_g,gauss);
    8596                 dhdt_input->GetInputValue(&dhdt_g,gauss);
    8597 
    8598                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8599                 GetNodalFunctions(basis,gauss);
    8600 
    8601                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
    8602         }
    8603 
    8604         /*Clean up and return*/
    8605         xDelete<IssmDouble>(basis);
    8606         delete gauss;
    8607         return pe;
    8608 }
    8609 /*}}}*/
    8610 /*FUNCTION Tria::CreatePVectorBalancethickness_DG {{{*/
    8611 ElementVector* Tria::CreatePVectorBalancethickness_DG(void){
    8612 
    8613         /*Intermediaries */
    8614         IssmDouble xyz_list[NUMVERTICES][3];
    8615         IssmDouble mb_g,ms_g,dhdt_g,Jdet;
    8616 
    8617         /*Fetch number of nodes and dof for this finite element*/
    8618         int numnodes = this->NumberofNodes();
    8619 
    8620         /*Initialize Element vector and other vectors*/
    8621         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    8622         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    8623 
    8624         /*Retrieve all inputs and parameters*/
    8625         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8626         Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
    8627         Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
    8628         Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);_assert_(dhdt_input);
    8629 
    8630         /* Start  looping on the number of gaussian points: */
    8631         GaussTria* gauss=new GaussTria(2);
    8632         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8633 
    8634                 gauss->GaussPoint(ig);
    8635 
    8636                 ms_input->GetInputValue(&ms_g,gauss);
    8637                 mb_input->GetInputValue(&mb_g,gauss);
    8638                 dhdt_input->GetInputValue(&dhdt_g,gauss);
    8639 
    8640                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8641                 GetNodalFunctions(basis,gauss);
    8642 
    8643                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*basis[i];
    8644         }
    8645 
    8646         /*Clean up and return*/
    8647         xDelete<IssmDouble>(basis);
    8648         delete gauss;
    8649         return pe;
    8650 }
    8651 /*}}}*/
    8652 /*FUNCTION Tria::CreatePVectorBalancevelocity{{{*/
    8653 ElementVector* Tria::CreatePVectorBalancevelocity(void){
    8654 
    8655         /*Intermediaries */
    8656         IssmDouble xyz_list[NUMVERTICES][3];
    8657         IssmDouble dhdt_g,mb_g,ms_g,Jdet;
    8658         IssmDouble h,gamma,thickness;
    8659         IssmDouble hnx,hny,dhnx[2],dhny[2];
    8660 
    8661         /*Fetch number of nodes and dof for this finite element*/
    8662         int numnodes = this->NumberofNodes();
    8663 
    8664         /*Initialize Element vector and other vectors*/
    8665         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    8666         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    8667         IssmDouble*    dbasis = xNew<IssmDouble>(numnodes*2);
    8668         IssmDouble*    H      = xNew<IssmDouble>(numnodes);
    8669         IssmDouble*    Nx     = xNew<IssmDouble>(numnodes);
    8670         IssmDouble*    Ny     = xNew<IssmDouble>(numnodes);
    8671 
    8672         /*Retrieve all inputs and parameters*/
    8673         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8674         Input* ms_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(ms_input);
    8675         Input* mb_input=inputs->GetInput(BasalforcingsMeltingRateEnum);   _assert_(mb_input);
    8676         Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum); _assert_(dhdt_input);
    8677         Input* H_input=inputs->GetInput(ThicknessEnum); _assert_(H_input);
    8678         h=sqrt(2.*this->GetArea());
    8679 
    8680         /*Get vector N for all nodes*/
    8681         GetInputListOnNodes(Nx,SurfaceSlopeXEnum);
    8682         GetInputListOnNodes(Ny,SurfaceSlopeYEnum);
    8683         GetInputListOnNodes(H,ThicknessEnum);
    8684         for(int i=0;i<numnodes;i++){
    8685                 IssmDouble norm=sqrt(Nx[i]*Nx[i]+Ny[i]*Ny[i]+1.e-10);
    8686                 Nx[i] = -H[i]*Nx[i]/norm;
    8687                 Ny[i] = -H[i]*Ny[i]/norm;
    8688         }
    8689 
    8690         /* Start  looping on the number of gaussian points: */
    8691         GaussTria* gauss=new GaussTria(2);
    8692         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8693 
    8694                 gauss->GaussPoint(ig);
    8695 
    8696                 ms_input->GetInputValue(&ms_g,gauss);
    8697                 mb_input->GetInputValue(&mb_g,gauss);
    8698                 dhdt_input->GetInputValue(&dhdt_g,gauss);
    8699                 H_input->GetInputValue(&thickness,gauss);
    8700                 if(thickness<50.) thickness=50.;
    8701 
    8702                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8703                 GetNodalFunctions(basis,gauss);
    8704                 GetNodalFunctionsDerivatives(dbasis,&xyz_list[0][0],gauss);
    8705 
    8706                 TriaRef::GetInputDerivativeValue(&dhnx[0],Nx,&xyz_list[0][0],gauss);
    8707                 TriaRef::GetInputDerivativeValue(&dhny[0],Ny,&xyz_list[0][0],gauss);
    8708                 TriaRef::GetInputValue(&hnx,Nx,gauss);
    8709                 TriaRef::GetInputValue(&hny,Ny,gauss);
    8710 
    8711                 gamma=h/(2.*thickness+1.e-10);
    8712 
    8713                 for(int i=0;i<numnodes;i++){
    8714                         pe->values[i]+=Jdet*gauss->weight*(ms_g-mb_g-dhdt_g)*(
    8715                                                 basis[i] + gamma*(basis[i]*(dhnx[0]+dhny[1])+hnx*dbasis[0*numnodes+i] + hny*dbasis[1*numnodes+i])
    8716                                                 );
    8717                 }
    8718         }
    8719 
    8720         /*Clean up and return*/
    8721         xDelete<IssmDouble>(basis);
    8722         xDelete<IssmDouble>(dbasis);
    8723         xDelete<IssmDouble>(H);
    8724         xDelete<IssmDouble>(Nx);
    8725         xDelete<IssmDouble>(Ny);
    8726         delete gauss;
    8727         return pe;
    8728 }
    8729 /*}}}*/
    8730 /*FUNCTION Tria::CreatePVectorSmoothedSlopeX{{{*/
    8731 ElementVector* Tria::CreatePVectorSmoothedSlopeX(void){
    8732 
    8733         /*Intermediaries */
    8734         IssmDouble xyz_list[NUMVERTICES][3];
    8735         IssmDouble Jdet;
    8736         IssmDouble thickness,slope[2];
    8737         IssmDouble taud_x,norms,normv,vx,vy;
    8738 
    8739         /*Fetch number of nodes and dof for this finite element*/
    8740         int numnodes = this->NumberofNodes();
    8741 
    8742         /*Initialize Element vector and other vectors*/
    8743         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    8744         IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    8745 
    8746         /*Retrieve all inputs and parameters*/
    8747         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8748         Input* H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
    8749         Input* surface_input = inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
    8750         Input* vx_input      = inputs->GetInput(VxEnum);
    8751         Input* vy_input      = inputs->GetInput(VyEnum);
    8752 
    8753         /* Start  looping on the number of gaussian points: */
    8754         GaussTria* gauss=new GaussTria(2);
    8755         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8756 
    8757                 gauss->GaussPoint(ig);
    8758 
    8759                 H_input->GetInputValue(&thickness,gauss);
    8760                 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    8761                 if(vx_input && vy_input){
    8762                         vx_input->GetInputValue(&vx,gauss);
    8763                         vy_input->GetInputValue(&vy,gauss);
    8764                         norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10);
    8765                         normv = sqrt(vx*vx + vy*vy);
    8766                         if(normv>15./(365.*24.*3600.)) slope[0] = -vx/normv*norms;
    8767                 }
    8768                 taud_x = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[0];
    8769 
    8770                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8771                 GetNodalFunctions(basis,gauss);
    8772 
    8773                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_x*basis[i];
    8774         }
    8775 
    8776         /*Clean up and return*/
    8777         xDelete<IssmDouble>(basis);
    8778         delete gauss;
    8779         return pe;
    8780 }
    8781 /*}}}*/
    8782 /*FUNCTION Tria::CreatePVectorSmoothedSlopeY{{{*/
    8783 ElementVector* Tria::CreatePVectorSmoothedSlopeY(void){
    8784 
    8785         /*Intermediaries */
    8786         IssmDouble xyz_list[NUMVERTICES][3];
    8787         IssmDouble Jdet;
    8788         IssmDouble thickness,slope[2];
    8789         IssmDouble taud_y,norms,normv,vx,vy;
    8790 
    8791         /*Fetch number of nodes and dof for this finite element*/
    8792         int numnodes = this->NumberofNodes();
    8793 
    8794         /*Initialize Element vector and other vectors*/
    8795         ElementVector* pe    = new ElementVector(nodes,numnodes,this->parameters);
    8796         IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    8797 
    8798         /*Retrieve all inputs and parameters*/
    8799         ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    8800         Input* H_input       = inputs->GetInput(ThicknessEnum); _assert_(H_input);
    8801         Input* surface_input = inputs->GetInput(SurfaceEnum);   _assert_(surface_input);
    8802         Input* vx_input      = inputs->GetInput(VxEnum);
    8803         Input* vy_input      = inputs->GetInput(VyEnum);
    8804 
    8805         /* Start  looping on the number of gaussian points: */
    8806         GaussTria* gauss=new GaussTria(2);
    8807         for(int ig=gauss->begin();ig<gauss->end();ig++){
    8808 
    8809                 gauss->GaussPoint(ig);
    8810 
    8811                 H_input->GetInputValue(&thickness,gauss);
    8812                 surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
    8813                 if(vx_input && vy_input){
    8814                         vx_input->GetInputValue(&vx,gauss);
    8815                         vy_input->GetInputValue(&vy,gauss);
    8816                         norms = sqrt(slope[0]*slope[0]+slope[1]*slope[1]+1.e-10);
    8817                         normv = sqrt(vx*vx + vy*vy);
    8818                         if(normv>15./(365.*24.*3600.)) slope[1] = -vy/normv*norms;
    8819                 }
    8820                 taud_y = matpar->GetRhoIce()*matpar->GetG()*thickness*slope[1];
    8821 
    8822                 GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
    8823                 GetNodalFunctions(basis,gauss);
    8824 
    8825                 for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*taud_y*basis[i];
    8826         }
    8827 
    8828         /*Clean up and return*/
    8829         xDelete<IssmDouble>(basis);
    8830         delete gauss;
    8831         return pe;
    88327132}
    88337133/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Tria.h

    r16956 r16982  
    222222                ElementMatrix* CreateMassMatrix(void);
    223223                ElementMatrix* CreateBasalMassMatrix(void);
    224                 ElementVector* CreatePVector(void);
    225                 ElementVector* CreatePVectorBalancethickness(void);
    226                 ElementVector* CreatePVectorBalancethickness_DG(void);
    227                 ElementVector* CreatePVectorBalancethickness_CG(void);
    228                 ElementVector* CreatePVectorBalancevelocity(void);
    229                 ElementVector* CreatePVectorSmoothedSlopeX(void);
    230                 ElementVector* CreatePVectorSmoothedSlopeY(void);
    231                 ElementVector* CreatePVectorMasstransport(void);
    232                 ElementVector* CreatePVectorMasstransport_CG(void);
    233                 ElementVector* CreatePVectorMasstransport_DG(void);
    234                 ElementVector* CreatePVectorFreeSurfaceTop(void);
    235                 ElementVector* CreatePVectorFreeSurfaceTop1D(void);
    236                 ElementVector* CreatePVectorFreeSurfaceBase(void);
    237                 ElementVector* CreatePVectorFreeSurfaceBase1D(void);
    238                 ElementVector* CreatePVectorL2Projection(void);
    239                 ElementVector* CreatePVectorL2ProjectionBase(void);
    240224                IssmDouble     EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};
    241225                IssmDouble     EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};
    … …  
    306290                ElementMatrix* CreateKMatrixStressbalanceFSViscous(void);
    307291                ElementMatrix* CreateKMatrixStressbalanceFSFriction(void);
    308                 ElementVector* CreatePVectorStressbalanceSSA(void);
    309                 ElementVector* CreatePVectorStressbalanceSSADrivingStress(void);
    310                 ElementVector* CreatePVectorStressbalanceSSAFront(void);
    311                 ElementVector* CreatePVectorStressbalanceSIA(void);
    312                 ElementVector* CreatePVectorStressbalanceFS(void);
    313                 ElementVector* CreatePVectorStressbalanceFSFront(void);
    314                 ElementVector* CreatePVectorStressbalanceFSViscous(void);
    315292                void           PVectorGLSstabilization(ElementVector* pe);
    316                 ElementVector* CreatePVectorStressbalanceFSShelf(void);
    317293                ElementMatrix* CreateJacobianStressbalanceSSA(void);
    318294                IssmDouble     GetYcoord(GaussTria* gauss);
    … …  
    322298                ElementMatrix* CreateKMatrixAdjointBalancethickness(void);
    323299                ElementMatrix* CreateKMatrixAdjointSSA(void);
    324                 ElementVector* CreatePVectorAdjointHoriz(void);
    325                 ElementVector* CreatePVectorAdjointBalancethickness(void);
    326300                #endif
    327301
    … …  
    339313                ElementMatrix* CreateKMatrixHydrologyDCEfficient(void);
    340314                ElementMatrix* CreateEPLDomainMassMatrix(void);
    341                 ElementVector* CreatePVectorHydrologyShreve(void);
    342                 ElementVector* CreatePVectorHydrologyDCInefficient(void);
    343                 ElementVector* CreatePVectorHydrologyDCEfficient(void);
    344                 ElementVector* CreatePVectorL2ProjectionEPL(void);
    345315                void           CreateHydrologyWaterVelocityInput(void);
    346316                void           GetHydrologyDCInefficientHmax(IssmDouble* ph_max,int index);
    … …  
    356326                ElementMatrix* CreateKMatrixDamageEvolution(void);
    357327                ElementMatrix* CreateKMatrixDamageEvolution_CG(void);
    358                 ElementVector* CreatePVectorDamageEvolution(void);
    359                 ElementVector* CreatePVectorDamageEvolution_CG(void);
    360328                void           DamageEvolutionF(IssmDouble* flist);
    361329                #endif

  • issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp

    r16819 r16982  
    4242                        element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    4343                        ElementMatrix* Ke = element->CreateKMatrix();
    44                         ElementVector* pe = element->CreatePVector();
     44                        //ElementVector* pe = element->CreatePVector();
     45                        ElementVector* pe = analysis->CreatePVector(element);
    4546                        element->ReduceMatrices(Ke,pe);
    4647                        if(Ke) Ke->AddToGlobal(Kff_temp,NULL);
    … …  
    7273        for (i=0;i<femmodel->elements->Size();i++){
    7374                element=dynamic_cast<Element*>(femmodel->elements->GetObjectByOffset(i));
    74                 //ElementVector* pe = analysis->CreatePVector(element);
     75                ElementMatrix* Ke = element->CreateKMatrix();
    7576                //ElementMatrix* Ke = analysis->CreateKMatrix(element);
    76                 ElementVector* pe = element->CreatePVector();
    77                 ElementMatrix* Ke = element->CreateKMatrix();
     77                //ElementVector* pe = element->CreatePVector();
     78                ElementVector* pe = analysis->CreatePVector(element);
    7879                element->ReduceMatrices(Ke,pe);
    7980                if(Ke) Ke->AddToGlobal(Kff,Kfs);
Note: See TracChangeset for help on using the changeset viewer.