Changeset 18502

Timestamp:

09/10/14 22:14:07 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added UpdateConstraints

Location:

issm/trunk-jpl/src/c

Files:

• analyses/SeaiceAnalysis.cpp (modified) (2 diffs)
• analyses/SeaiceAnalysis.h (modified) (1 diff)
• cores/seaice_core.cpp (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.cpp

@@ -321 +321 @@
 }/*}}}*/
 void SeaiceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
-        /*Default, do nothing*/
-        return;
+
+        /*Intermediaries*/
+        Vector<IssmDouble>* mask        = NULL;
+        IssmDouble*         serial_mask = NULL;
+
+        /*Step 1: update mask of active nodes*/
+        mask=new Vector<IssmDouble>(femmodel->nodes->NumberOfNodes(SeaiceAnalysisEnum));
+        for (int i=0;i<femmodel->elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(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;in<numnodes;in++) mask->SetValue(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;i<femmodel->nodes->Size();i++){
+                Node* node=dynamic_cast<Node*>(femmodel->nodes->GetObjectByOffset(i));
+                if(node->InAnalysis(SeaiceAnalysisEnum)){
+                        if(serial_mask[node->Sid()]==1.){
+                                node->Activate();
+                                counter++;
+                        }
+                        else{
+                                node->Deactivate();
+                        }
+                }
+        }
+        xDelete<IssmDouble>(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 +463 @@
         xDelete<IssmDouble>(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;i<femmodel->elements->Size();i++){
+                Element* element=dynamic_cast<Element*>(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);
+        }
+}/*}}}*/

issm/trunk-jpl/src/c/analyses/SeaiceAnalysis.h

@@ -32 +32 @@
 
                 /*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);

issm/trunk-jpl/src/c/cores/seaice_core.cpp

@@ -18 +18 @@
         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);