Index: /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp (revision 18504) +++ /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp (revision 18505) @@ -141,5 +141,5 @@ for(int i=0;i<3;i++){ for(int j=0;j<3;j++){ - C[i][j] = dt*time_ratio*gauss->weight*Jdet*C[i][j]; + C[i][j] = dt*thickness*time_ratio*gauss->weight*Jdet*C[i][j]; } } @@ -320,6 +320,5 @@ /*Intermediaries*/ - int* doflist=NULL; - IssmDouble* xyz_list=NULL; + int* doflist=NULL; /*Fetch number of nodes and dof for this finite element*/ @@ -358,5 +357,4 @@ xDelete(vx); xDelete(values); - xDelete(xyz_list); xDelete(doflist); }/*}}}*/ @@ -416,5 +414,5 @@ /*Intermediaries*/ IssmDouble D[3][3]; - IssmDouble damage,thickness,concentration; + IssmDouble damage,concentration; /*Get Material parameters*/ @@ -425,8 +423,6 @@ /*Get damage input at this location*/ Input* damage_input = element->GetInput(DamageEnum); _assert_(damage_input); - Input* thickness_input = element->GetInput(SeaiceThicknessEnum); _assert_(thickness_input); Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input); damage_input->GetInputValue(&damage,gauss); - thickness_input->GetInputValue(&thickness,gauss); concentration_input->GetInputValue(&concentration,gauss); @@ -448,5 +444,5 @@ for(int i=0;i<3;i++){ for(int j=0;j<3;j++){ - C[i*3+j] = thickness*E*D[i][j]; + C[i*3+j] = E*D[i][j]; } } @@ -623,2 +619,150 @@ } }/*}}}*/ +void SeaiceAnalysis::PostProcess(FemModel* femmodel){/*{{{*/ + + /*Intermediaties*/ + IssmDouble *xyz_list = NULL; + IssmDouble *xyz_list_old = NULL; + IssmDouble *xyz_list_new = NULL; + IssmDouble dt,meshx,meshy,vx,vy,eps_xx,eps_yy,eps_zz,area_old,area_new; + IssmDouble time_relaxation_stress,time_relaxation_damage,damage,concentration,thickness; + IssmDouble sigma_dot_xx,sigma_dot_yy,sigma_dot_xy; + IssmDouble sigma_pred_xx,sigma_pred_yy,sigma_pred_xy; + IssmDouble sigma_xx,sigma_yy,sigma_xy; + IssmDouble min_c,min_h,max_h; + IssmDouble epsilon[3]; /*eps_xx,eps_yy,eps_xy*/ + IssmDouble C[3][3]; + + /*Fetch the parameters once for all*/ + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&min_c,SeaiceMinConcentrationEnum); + femmodel->parameters->FindParam(&min_h,SeaiceMinThicknessEnum); + femmodel->parameters->FindParam(&max_h,SeaiceMaxThicknessEnum); + + /*Loop over the elements of this partition and update accordingly*/ + for(int i=0;ielements->Size();i++){ + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + + /*Get some inputs needed for the update*/ + element->GetVerticesCoordinates(&xyz_list); + time_relaxation_stress = element->GetMaterialParameter(MaterialsTimeRelaxationStressEnum); + time_relaxation_damage = element->GetMaterialParameter(MaterialsTimeRelaxationDamageEnum); + Input* meshx_input = element->GetInput(MeshXEnum); _assert_(meshx_input); + Input* meshy_input = element->GetInput(MeshYEnum); _assert_(meshy_input); + Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); + + /*Preallocate future inputs*/ + int numvertices=element->GetNumberOfVertices(); + IssmDouble* meshx_new_list = xNew(numvertices); + IssmDouble* meshy_new_list = xNew(numvertices); + xyz_list_old = xNew(numvertices*3); + xyz_list_new = xNew(numvertices*3); + + /*1. update vertex positions (using a Gauss object for convenience)*/ + Gauss* gauss=element->NewGauss(); + for (int iv=0;ivGaussVertex(iv); + + meshx_input->GetInputValue(&meshx,gauss); + meshy_input->GetInputValue(&meshy,gauss); + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + + meshx_new_list[iv] = meshx + vx*dt; + meshy_new_list[iv] = meshy + vy*dt; + + xyz_list_old[iv*3+0] = meshx; + xyz_list_old[iv*3+1] = meshy; + xyz_list_old[iv*3+2] = 0.; + + xyz_list_new[iv*3+0] = meshx_new_list[iv]; + xyz_list_new[iv*3+1] = meshy_new_list[iv]; + xyz_list_new[iv*3+2] = 0.; + } + element->AddInput(MeshXEnum,meshx_new_list,P1Enum); + element->AddInput(MeshYEnum,meshy_new_list,P1Enum); + + /*Now we are going to use a point in the center of the element, because we are P0 for the stress*/ + gauss->GaussNode(P0Enum,0); + + /*Calculate sigma_dot*/ + sigma_dot_xx=sigma_dot_yy=sigma_dot_xy=0.; + this->CreateCTensor(&C[0][0],element,gauss); + element->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input); + for(int j=0;j<3;j++){ + sigma_dot_xx += C[0][j]*epsilon[0]; + sigma_dot_yy += C[1][j]*epsilon[1]; + sigma_dot_xy += C[2][j]*epsilon[2]; + } + + /*Get predicted stress state*/ + Input* sigma_xx_input = element->GetInput(StressTensorPredictorxxEnum); _assert_(sigma_xx_input); + Input* sigma_yy_input = element->GetInput(StressTensorPredictoryyEnum); _assert_(sigma_yy_input); + Input* sigma_xy_input = element->GetInput(StressTensorPredictorxyEnum); _assert_(sigma_xy_input); + sigma_xx_input->GetInputAverage(&sigma_xx); + sigma_yy_input->GetInputAverage(&sigma_yy); + sigma_xy_input->GetInputAverage(&sigma_xy); + + /*Calculate new stress and push to element*/ + sigma_xx = sigma_xx + dt*(sigma_dot_xx - sigma_xx/time_relaxation_stress); + sigma_yy = sigma_yy + dt*(sigma_dot_yy - sigma_yy/time_relaxation_stress); + sigma_xy = sigma_xy + dt*(sigma_dot_xy - sigma_xy/time_relaxation_stress); + element->AddInput(StressTensorxxEnum,&sigma_xx,P0Enum); + element->AddInput(StressTensoryyEnum,&sigma_yy,P0Enum); + element->AddInput(StressTensorxyEnum,&sigma_xy,P0Enum); + + /*Update Damage According*/ + Input* damage_input = element->GetInput(DamageEnum); _assert_(damage_input); + damage_input->GetInputAverage(&damage); + damage = damage*(1.-dt/time_relaxation_damage); + element->AddInput(DamageEnum,&damage,P0Enum); + + /*Prepare new predictor*/ + this->CreateCTensor(&C[0][0],element,gauss); /*C is updated now that there is a new damage in inputs*/ + for(int j=0;j<3;j++){ + sigma_pred_xx += C[0][j]*epsilon[0]; + sigma_pred_yy += C[1][j]*epsilon[1]; + sigma_pred_xy += C[2][j]*epsilon[2]; + } + element->AddInput(StressTensorPredictorxxEnum,&sigma_pred_xx,P0Enum); + element->AddInput(StressTensorPredictoryyEnum,&sigma_pred_yy,P0Enum); + element->AddInput(StressTensorPredictorxyEnum,&sigma_pred_xy,P0Enum); + + /*Calculate Old and new area (FIXME: for now we assume trianlges...)*/ + element->JacobianDeterminant(&area_old,xyz_list_old,gauss); + element->JacobianDeterminant(&area_new,xyz_list_new,gauss); + area_old = area_old/SQRT3; + area_new = area_new/SQRT3; + + /*Update ice thickness and concentration using element distortion*/ + Input* thickness_input = element->GetInput(SeaiceThicknessEnum); _assert_(thickness_input); + Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input); + thickness_input->GetInputAverage(&thickness); + concentration_input->GetInputAverage(&concentration); + IssmDouble ice_area = concentration*area_old; + IssmDouble ice_volume = thickness*area_old; + if(concentration>min_c){ + concentration = ice_area /area_new; + thickness = ice_volume/area_new; + + /* lower bounds */ + concentration = ((concentration>min_c)?(concentration):(min_c)); + thickness = ((thickness >min_h)?(thickness ):(min_h)); + + /* upper bounds (only for the concentration) */ + concentration = ((concentration<1. )?(concentration):(1)); + thickness = ((thickness AddInput(SeaiceConcentrationEnum,&concentration,P0Enum); + element->AddInput(SeaiceThicknessEnum,&thickness,P0Enum); + + /*Clean up*/ + xDelete(meshx_new_list); + xDelete(meshy_new_list); + xDelete(xyz_list); + xDelete(xyz_list_old); + xDelete(xyz_list_new); + } + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h (revision 18504) +++ /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h (revision 18505) @@ -36,4 +36,5 @@ void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void GetM(IssmDouble* M,Element* element,IssmDouble* xyz_list,Gauss* gauss); + void PostProcess(FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/cores/seaice_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/seaice_core.cpp (revision 18504) +++ /issm/trunk-jpl/src/c/cores/seaice_core.cpp (revision 18505) @@ -13,5 +13,6 @@ /*parameters: */ - bool save_results; + bool save_results; + SeaiceAnalysis analysis = SeaiceAnalysis(); /*recover parameters: */ @@ -20,12 +21,16 @@ /*Launch solution sequence for the only analysis we are interested in*/ femmodel->SetCurrentConfiguration(SeaiceAnalysisEnum); - SeaiceAnalysis* analysis = new SeaiceAnalysis(); - analysis->UpdateDamageAndStress(femmodel); + analysis.UpdateDamageAndStress(femmodel); solutionsequence_linear(femmodel); + analysis.PostProcess(femmodel); if(save_results){ if(VerboseSolution()) _printf0_(" saving results\n"); - int outputs[6] = {VxEnum,VyEnum,VelEnum,SeaiceConcentrationEnum,SeaiceThicknessEnum,DamageEnum}; - femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],6); + int outputs[14] = {VxEnum,VyEnum,VelEnum, + SeaiceConcentrationEnum,SeaiceThicknessEnum,DamageEnum, + StressTensorPredictorxxEnum,StressTensorPredictoryyEnum,StressTensorPredictorxyEnum, + StressTensorxxEnum,StressTensoryyEnum,StressTensorxyEnum, + MeshXEnum,MeshYEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],14); }