Index: /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp (revision 18501) +++ /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp (revision 18502) @@ -321,6 +321,47 @@ }/*}}}*/ void SeaiceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ - /*Default, do nothing*/ - return; + + /*Intermediaries*/ + Vector* mask = NULL; + IssmDouble* serial_mask = NULL; + + /*Step 1: update mask of active nodes*/ + mask=new Vector(femmodel->nodes->NumberOfNodes(SeaiceAnalysisEnum)); + for (int i=0;ielements->Size();i++){ + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + + /*Get current concentration of element and decide whether it is an active element*/ + int numnodes = element->GetNumberOfNodes(); + Input* concentration_input = element->GetInput(SeaiceConcentrationEnum); _assert_(concentration_input); + if(concentration_input->Max()>0.){ + for(int in=0;inSetValue(element->nodes[in]->Sid(),1.,INS_VAL); + } + } + + /*Assemble and serialize*/ + mask->Assemble(); + serial_mask=mask->ToMPISerial(); + delete mask; + + /*Update node activation accordingly*/ + int counter =0; + for(int i=0;inodes->Size();i++){ + Node* node=dynamic_cast(femmodel->nodes->GetObjectByOffset(i)); + if(node->InAnalysis(SeaiceAnalysisEnum)){ + if(serial_mask[node->Sid()]==1.){ + node->Activate(); + counter++; + } + else{ + node->Deactivate(); + } + } + } + xDelete(serial_mask); + + /*Display number of active nodes*/ + int sum_counter; + ISSM_MPI_Reduce(&counter,&sum_counter,1,ISSM_MPI_INT,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + if(VerboseSolution()) _printf0_(" Number of active nodes: "<< sum_counter <<"\n"); }/*}}}*/ @@ -422,2 +463,121 @@ xDelete(basis); }/*}}}*/ +void SeaiceAnalysis::UpdateDamageAndStress(FemModel* femmodel){/*{{{*/ + /* The damage variable is updated as a function of the actual elastic deformation + * In both cases, a Coulombic enveloppe is used, define by the cohesion C, tan(phi) and tract_coef. + * In both cases, a maximal compressive strength is fixed at compr_max + * The enveloppe is defined in N/m^2. + * The coeficients of the internal stress are then multiplied by the ice thickness to be used in the vertical integrated momentiom equation. + */ + + /* Mohr-Coulomb criterion calculation + * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + * + * sigma_s + * Mohr.Coulomb branch | + * \ | + * * | | + * * | cohesion (C=cfix+calea) + * | * | / + * | * | / tract + * | * |/ / + * -comp_max | * / + * \ | | * / + * \| 0| | * + * -------------------------*------------ sigma_n + * | | | * + * | | * + * | * + * | * | + * | * | + * | * | + * * | + * * | | + * * | + * | + * | + * + * %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + */ + + /*Intermediaties*/ + IssmDouble sigma_xx,sigma_yy,sigma_xy,sigma_s,sigma_n,sigma_target; + IssmDouble compression_coef,traction_coef,cohesion,tan_phi; + IssmDouble traction,compression_max; + IssmDouble damage_test,damage_new,damage,tmp; + + /*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 Mohr-Coulomb parameters*/ + compression_coef = element->GetMaterialParameter(MaterialsCompressionCoefEnum); + traction_coef = element->GetMaterialParameter(MaterialsTractionCoefEnum); + cohesion = element->GetMaterialParameter(MaterialsCohesionEnum); //C + tan_phi = element->GetMaterialParameter(MaterialsInternalFrictionCoefEnum); //mu = tan(phi) + + /*Get current stress state*/ + Input* sigma_xx_input = element->GetInput(StressTensorxxEnum); + Input* sigma_yy_input = element->GetInput(StressTensoryyEnum); + Input* sigma_xy_input = element->GetInput(StressTensorxyEnum); + if(sigma_xx_input) sigma_xx_input->GetInputAverage(&sigma_xx); + else sigma_xx = 0.; + if(sigma_yy_input) sigma_yy_input->GetInputAverage(&sigma_yy); + else sigma_yy = 0.; + if(sigma_xy_input) sigma_xy_input->GetInputAverage(&sigma_xy); + else sigma_xy = 0.; + + /* Compute the invariants of the elastic deformation and instantaneous deformation rate */ + sigma_s=sqrt(pow((sigma_xx-sigma_yy)/2.,2)+pow(sigma_xy,2)); + sigma_n=(sigma_xx+sigma_yy)/2.; + + /* same maximal tensile strength */ + traction=traction_coef*cohesion/tan_phi; + compression_max=compression_coef*cohesion; + + /* estimate the internal constraints using the current elastic deformation */ + Input* damage_input = element->GetInput(MaterialsDamageEnum); _assert_(damage_input); + Input* damagenew_input = element->GetInput(MaterialsDamageEnum); _assert_(damagenew_input);//FIXME + damage_input->GetInputAverage(&damage); + damagenew_input->GetInputAverage(&damage_new); + damage_test = damage; + if(sigma_n>traction || sigma_n<-compression_max){ + if(sigma_n>traction){ + sigma_target=traction; + } + else{ + sigma_target=-compression_max; + } + + tmp=1.-sigma_target/sigma_n*(1.-damage); + if(tmp<1. && tmp>damage_new){ + damage_test=((damage_test>tmp)?(damage_test):(tmp)); /* max(damage_test,tmp); */ + } + } + if(sigma_s>cohesion-sigma_n*tan_phi){ + tmp=1.0-cohesion/(sigma_s+sigma_n*tan_phi)*(1-damage); + if(tmp<1. && tmp>damage_new){ + damage_test=((damage_test>tmp)?(damage_test):(tmp)); /*max(damage_test,tmp); */ + } + } + + /* The damage variable is changed */ + damage_new=damage_test; + element->AddInput(MaterialsDamageEnum,&damage_test,P0Enum);//FIXME + + /* Recompute the internal stress*/ + if(damage<1.){ + sigma_xx = (1.-damage_new)/(1.-damage)*sigma_xx; + sigma_yy = (1.-damage_new)/(1.-damage)*sigma_yy; + sigma_xy = (1.-damage_new)/(1.-damage)*sigma_xy; + } + else{ + sigma_xx = 0.; + sigma_yy = 0.; + sigma_xy = 0.; + } + element->AddInput(StressTensorxxEnum,&sigma_xx,P0Enum); + element->AddInput(StressTensoryyEnum,&sigma_yy,P0Enum); + element->AddInput(StressTensorxyEnum,&sigma_xy,P0Enum); + } +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h (revision 18501) +++ /issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h (revision 18502) @@ -32,4 +32,5 @@ /*Sea ice specifics*/ + void UpdateDamageAndStress(FemModel* femmodel); void CreateCTensor(IssmDouble* C,Element* element,Gauss* gauss); void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); Index: /issm/trunk-jpl/src/c/cores/seaice_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/seaice_core.cpp (revision 18501) +++ /issm/trunk-jpl/src/c/cores/seaice_core.cpp (revision 18502) @@ -18,5 +18,5 @@ femmodel->parameters->FindParam(&save_results,SaveResultsEnum); - if(VerboseSolution()) _printf0_("call computational core:\n"); + /*Launch solution sequence for the only analysis we are interested in*/ femmodel->SetCurrentConfiguration(SeaiceAnalysisEnum); solutionsequence_linear(femmodel);