Index: /issm/trunk-jpl/src/c/analyses/AgeAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AgeAnalysis.cpp (revision 27571) +++ /issm/trunk-jpl/src/c/analyses/AgeAnalysis.cpp (revision 27572) @@ -92,68 +92,4 @@ }/*}}}*/ ElementMatrix* AgeAnalysis::CreateKMatrix(Element* element){/*{{{*/ - - _error_("STOP"); - /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - - /*compute all stiffness matrices for this element*/ - ElementMatrix* Ke1=CreateKMatrixVolume(element); - ElementMatrix* Ke2=CreateKMatrixShelf(element); - ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); - - /*clean-up and return*/ - delete Ke1; - delete Ke2; - return Ke; -}/*}}}*/ -ElementMatrix* AgeAnalysis::CreateKMatrixShelf(Element* element){/*{{{*/ - - /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - - /*Initialize Element matrix and return if necessary*/ - if(!element->IsOnBase() || !element->IsAllFloating()) return NULL; - - IssmDouble dt,Jdet,D; - IssmDouble *xyz_list_base = NULL; - - /*Get basal element*/ - if(!element->IsOnBase() || !element->IsAllFloating()) return NULL; - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Initialize vectors*/ - ElementMatrix* Ke = element->NewElementMatrix(); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - element->FindParam(&dt,TimesteppingTimeStepEnum); - IssmDouble gravity = element->FindParam(ConstantsGEnum); - IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum); - IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); - IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum); - IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGaussBase(4); - while(gauss->next()){ - - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctions(basis,gauss); - - D=gauss->weight*Jdet*rho_water*mixed_layer_capacity/(heatcapacity*rho_ice); - if(reCast(dt)) D=dt*D; - for(int i=0;ivalues[i*numnodes+j] += D*basis[i]*basis[j]; - } - - /*Clean up and return*/ - delete gauss; - xDelete(basis); - xDelete(xyz_list_base); - return Ke; -}/*}}}*/ -ElementMatrix* AgeAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/ /* Check if ice in element */ @@ -303,135 +239,7 @@ if(!element->IsIceInElement()) return NULL; - /*compute all load vectors for this element*/ - ElementVector* pe1=CreatePVectorVolume(element); - ElementVector* pe2=CreatePVectorSheet(element); - ElementVector* pe3=CreatePVectorShelf(element); - ElementVector* pe =new ElementVector(pe1,pe2,pe3); - - /*clean-up and return*/ - delete pe1; - delete pe2; - delete pe3; - return pe; -}/*}}}*/ -ElementVector* AgeAnalysis::CreatePVectorSheet(Element* element){/*{{{*/ - - /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - - /* Geothermal flux on ice sheet base and basal friction */ - if(!element->IsOnBase() || element->IsAllFloating()) return NULL; - - IssmDouble dt,Jdet,geothermalflux,vx,vy,vz; - IssmDouble alpha2,scalar,basalfriction,heatflux; - IssmDouble *xyz_list_base = NULL; - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Initialize vectors*/ - ElementVector* pe = element->NewElementVector(); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - element->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input = element->GetInput(VzEnum); _assert_(vz_input); - Input* geothermalflux_input = element->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input); - IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); - IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum); - - /*Build friction element, needed later: */ - Friction* friction=new Friction(element,3); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss = element->NewGaussBase(4); - while(gauss->next()){ - - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctions(basis,gauss); - - geothermalflux_input->GetInputValue(&geothermalflux,gauss); - friction->GetAlpha2(&alpha2,gauss); - vx_input->GetInputValue(&vx,gauss); - vy_input->GetInputValue(&vy,gauss); - vz_input->GetInputValue(&vz,gauss); - vz = 0.;//FIXME - basalfriction = alpha2*(vx*vx + vy*vy + vz*vz); - heatflux = (basalfriction+geothermalflux)/(rho_ice*heatcapacity); - - scalar = gauss->weight*Jdet*heatflux; - if(dt!=0.) scalar=dt*scalar; - - for(int i=0;ivalues[i]+=scalar*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - delete friction; - xDelete(basis); - xDelete(xyz_list_base); - return pe; -}/*}}}*/ -ElementVector* AgeAnalysis::CreatePVectorShelf(Element* element){/*{{{*/ - - /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - - IssmDouble t_pmp,dt,Jdet,scalar_ocean,pressure; - IssmDouble *xyz_list_base = NULL; - - /*Get basal element*/ - if(!element->IsOnBase() || !element->IsAllFloating()) return NULL; - - /*Fetch number of nodes for this finite element*/ - int numnodes = element->GetNumberOfNodes(); - - /*Initialize vectors*/ - ElementVector* pe = element->NewElementVector(); - IssmDouble* basis = xNew(numnodes); - - /*Retrieve all inputs and parameters*/ - element->GetVerticesCoordinatesBase(&xyz_list_base); - element->FindParam(&dt,TimesteppingTimeStepEnum); - Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input); - IssmDouble gravity = element->FindParam(ConstantsGEnum); - IssmDouble rho_water = element->FindParam(MaterialsRhoSeawaterEnum); - IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); - IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum); - IssmDouble mixed_layer_capacity= element->FindParam(MaterialsMixedLayerCapacityEnum); - - /* Start looping on the number of gaussian points: */ - Gauss* gauss=element->NewGaussBase(4); - while(gauss->next()){ - - element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); - element->NodalFunctions(basis,gauss); - - pressure_input->GetInputValue(&pressure,gauss); - t_pmp=element->TMeltingPoint(pressure); - - scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*(t_pmp)/(heatcapacity*rho_ice); - if(reCast(dt)) scalar_ocean=dt*scalar_ocean; - - for(int i=0;ivalues[i]+=scalar_ocean*basis[i]; - } - - /*Clean up and return*/ - delete gauss; - xDelete(basis); - xDelete(xyz_list_base); - return pe; -}/*}}}*/ -ElementVector* AgeAnalysis::CreatePVectorVolume(Element* element){/*{{{*/ - - /* Check if ice in element */ - if(!element->IsIceInElement()) return NULL; - /*Intermediaries*/ int stabilization; - IssmDouble Jdet,phi,dt; + IssmDouble Jdet,dt; IssmDouble temperature; IssmDouble tau_parameter,diameter,hx,hy,hz; @@ -451,8 +259,4 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); - IssmDouble rho_ice = element->FindParam(MaterialsRhoIceEnum); - IssmDouble heatcapacity = element->FindParam(MaterialsHeatcapacityEnum); - IssmDouble thermalconductivity = 1.; - IssmDouble kappa = thermalconductivity/(rho_ice*heatcapacity); element->FindParam(&dt,TimesteppingTimeStepEnum); element->FindParam(&stabilization,AgeStabilizationEnum); @@ -469,7 +273,6 @@ element->JacobianDeterminant(&Jdet,xyz_list,gauss); element->NodalFunctions(basis,gauss); - element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input); - - scalar_def=phi/(rho_ice*heatcapacity)*Jdet*gauss->weight; + + scalar_def=1.*Jdet*gauss->weight; if(reCast(dt)) scalar_def=scalar_def*dt; @@ -490,5 +293,5 @@ vz_input->GetInputValue(&w,gauss); - tau_parameter=element->StabilizationParameter(u,v,w,diameter,kappa); + tau_parameter=element->StabilizationParameter(u,v,w,diameter,1.e-15); //assume very small conductivity to get tau for(int i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0*numnodes+i]+v*dbasis[1*numnodes+i]+w*dbasis[2*numnodes+i]); @@ -504,5 +307,5 @@ vy_input->GetInputValue(&v,gauss); vz_input->GetInputValue(&w,gauss); - element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,kappa); + element->StabilizationParameterAnisotropic(&tau_parameter_anisotropic[0],u,v,w,hx,hy,hz,1.e-15); //assume very small conductivity to get tau tau_parameter_hor=tau_parameter_anisotropic[0]; tau_parameter_ver=tau_parameter_anisotropic[1]; @@ -518,5 +321,4 @@ delete gauss; return pe; - }/*}}}*/ void AgeAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ @@ -533,5 +335,39 @@ void AgeAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ SetActiveNodesLSMx(femmodel); - _error_("Should also automatically constrain surface/basal nodes where we have inflow"); -}/*}}}*/ + + /*Constrain all nodes that are grounded and unconstrain the ones that float*/ + for(Object* & object : femmodel->elements->objects){ + Element *element = xDynamicCast(object); + + if(element->IsOnSurface()){ + element = element->SpawnTopElement(); + int numnodes = element->GetNumberOfNodes(); + IssmDouble *mask = xNew(numnodes); + IssmDouble *bed = xNew(numnodes); + IssmDouble *ls_active = xNew(numnodes); + + element->GetInputListOnNodes(&mask[0],MaskOceanLevelsetEnum); + element->GetInputListOnNodes(&bed[0],BaseEnum); + element->GetInputListOnNodes(&ls_active[0],IceMaskNodeActivationEnum); + + for(int in=0;inGetNode(in); + if(mask[in]<0. && ls_active[in]==1.){ + node->Activate(); + } + else{ + node->Deactivate(); + node->ApplyConstraint(0,bed[in]); + } + } + xDelete(mask); + xDelete(bed); + xDelete(ls_active); + } + else if(element->IsOnBase()){ + element = element->SpawnBasalElement(); + _error_("not supported"); + } + } +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/AgeAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/AgeAnalysis.h (revision 27571) +++ /issm/trunk-jpl/src/c/analyses/AgeAnalysis.h (revision 27572) @@ -26,10 +26,5 @@ ElementMatrix* CreateJacobianMatrix(Element* element); ElementMatrix* CreateKMatrix(Element* element); - ElementMatrix* CreateKMatrixShelf(Element* element); - ElementMatrix* CreateKMatrixVolume(Element* element); ElementVector* CreatePVector(Element* element); - ElementVector* CreatePVectorSheet(Element* element); - ElementVector* CreatePVectorShelf(Element* element); - ElementVector* CreatePVectorVolume(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); void GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index);