Changeset 16848 for issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

Timestamp:

11/21/13 09:21:32 (11 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with Damage evolution

File:

• issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (modified) (2 diffs)

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

@@ -99 +99 @@
 }/*}}}*/
 ElementVector* DamageEvolutionAnalysis::CreatePVector(Element* element){/*{{{*/
-_error_("not implemented yet");
+
+        /*Intermediaries*/
+        int      meshtype;
+        Element* basalelement;
+
+        /*Get basal element*/
+        element->FindParam(&meshtype,MeshTypeEnum);
+        switch(meshtype){
+                case Mesh2DhorizontalEnum:
+                        basalelement = element;
+                        break;
+                case Mesh3DEnum:
+                        if(!element->IsOnBed()) return NULL;
+                        basalelement = element->SpawnBasalElement();
+                        break;
+                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+        }
+
+        /*Intermediaries */
+        IssmDouble  Jdet,dt;
+        IssmDouble  f,damage;
+        IssmDouble* xyz_list = NULL;
+
+        /*Fetch number of nodes and dof for this finite element*/
+        int numnodes = element->GetNumberOfNodes();
+
+        /*Initialize Element vector and other vectors*/
+        ElementVector* pe    = element->NewElementVector();
+        IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+        /*Retrieve all inputs and parameters*/
+        element->GetVerticesCoordinates(&xyz_list);
+        element->FindParam(&dt,TimesteppingTimeStepEnum);
+        this->CreateDamageFInput(element);
+        Input* damaged_input = element->GetMaterialInput(DamageDbarEnum); _assert_(damaged_input);
+        Input* damagef_input = element->GetMaterialInput(DamageFEnum);    _assert_(damagef_input);
+
+        /* Start  looping on the number of gaussian points: */
+        Gauss* gauss=element->NewGauss(2);
+        for(int ig=gauss->begin();ig<gauss->end();ig++){
+                gauss->GaussPoint(ig);
+
+                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+                element->NodalFunctions(basis,gauss);
+
+                damaged_input->GetInputValue(&damage,gauss);
+                damagef_input->GetInputValue(&f,gauss);
+
+                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*(damage+dt*f)*basis[i];
+        }
+
+        /*Clean up and return*/
+        xDelete<IssmDouble>(xyz_list);
+        xDelete<IssmDouble>(basis);
+        if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+        delete gauss;
+        return pe;
 }/*}}}*/
 void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
@@ -136 +192 @@
         xDelete<int>(doflist);
 }/*}}}*/
+
+/*Intermediaries*/
+void DamageEvolutionAnalysis::CreateDamageFInput(Element* element){/*{{{*/
+
+        /*Intermediaries */
+        IssmDouble c1,c2,c3,healing,stress_threshold;
+        IssmDouble s_xx,s_xy,s_yy;
+        IssmDouble J2s,Xis,Psi,PosPsi,NegPsi;
+        IssmDouble damage,sigma_xx,sigma_xy,sigma_yy;
+
+        /*Fetch number of vertices and allocate output*/
+        int numvertices = element->GetNumberOfVertices();
+        IssmDouble* f   = xNew<IssmDouble>(numvertices);
+
+        /*retrieve parameters:*/
+        element->FindParam(&c1,DamageC1Enum);
+        element->FindParam(&c2,DamageC2Enum);
+        element->FindParam(&c3,DamageC3Enum);
+        element->FindParam(&healing,DamageHealingEnum);
+        element->FindParam(&stress_threshold,DamageStressThresholdEnum);
+
+        /*Compute stress tensor: */
+        element->ComputeStressTensor();
+
+        /*retrieve what we need: */
+        Input* sigma_xx_input  = element->GetInput(StressTensorxxEnum);     _assert_(sigma_xx_input);
+        Input* sigma_xy_input  = element->GetInput(StressTensorxyEnum);     _assert_(sigma_xy_input);
+        Input* sigma_yy_input  = element->GetInput(StressTensoryyEnum);     _assert_(sigma_yy_input);
+        Input* damage_input    = element->GetMaterialInput(DamageDbarEnum); _assert_(damage_input);
+
+        /*Damage evolution z mapping: */
+        Gauss* gauss=element->NewGauss();
+        for (int iv=0;iv<numvertices;iv++){
+                gauss->GaussVertex(iv);
+                
+                damage_input->GetInputValue(&damage,gauss);
+                sigma_xx_input->GetInputValue(&sigma_xx,gauss);
+                sigma_xy_input->GetInputValue(&sigma_xy,gauss);
+                sigma_yy_input->GetInputValue(&sigma_yy,gauss);
+
+                s_xx=sigma_xx/(1.-damage);
+                s_xy=sigma_xy/(1.-damage);
+                s_yy=sigma_yy/(1.-damage);
+
+                J2s=1./sqrt(2.)*sqrt(s_xx*s_xx + s_yy*s_yy + s_xy*s_xy);
+                Xis=sqrt(3.0)*J2s;
+                Psi=Xis-stress_threshold;
+                PosPsi=max(Psi,0.);
+                NegPsi=max(-Psi,0.);
+
+                f[iv]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1.-damage),-c3);
+        }
+
+        /*Add input*/
+        element->AddMaterialInput(DamageFEnum,f,P1Enum);
+        
+        /*Clean up and return*/
+        delete gauss;
+}/*}}}*/