source: issm/oecreview/Archive/19101-20495/ISSM-19387-19388.diff@ 20498

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

@@ -129 +129 @@
         IssmDouble epsf,stress_threshold,eps0;
         IssmDouble damage,B,n,epseff;
         IssmDouble eps_xx,eps_yy,eps_xy,eps1,eps2,epstmp;
-        int domaintype,damagelaw;
+        int domaintype;
 
         /*Fetch number of vertices and allocate output*/
         int numnodes = element->GetNumberOfNodes();
@@ -139 +139 @@
         element->FindParam(&epsf,DamageC1Enum);
         element->FindParam(&stress_threshold,DamageStressThresholdEnum);
         element->FindParam(&domaintype,DomainTypeEnum);
-        element->FindParam(&damagelaw,DamageLawEnum);
 
         /*Compute stress tensor: */
         element->ComputeStrainRate();
@@ -201 +200 @@
         IssmDouble s_xx,s_xy,s_xz,s_yy,s_yz,s_zz,s1,s2,s3,stmp;
         IssmDouble J2s,Chi,Psi,PosPsi,NegPsi;
         IssmDouble damage,tau_xx,tau_xy,tau_xz,tau_yy,tau_yz,tau_zz,stressMaxPrincipal;
-        int equivstress,domaintype,damagelaw,dim;
+        int equivstress,domaintype,dim;
 
         /*Fetch number of vertices and allocate output*/
         int numnodes = element->GetNumberOfNodes();
@@ -214 +213 @@
         element->FindParam(&healing,DamageHealingEnum);
         element->FindParam(&stress_threshold,DamageStressThresholdEnum);
         element->FindParam(&domaintype,DomainTypeEnum);
-        element->FindParam(&damagelaw,DamageLawEnum);
 
         /*Get problem dimension*/
         switch(domaintype){
@@ -290 +288 @@
                         }
                         Psi=Chi-stress_threshold;
                         NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
-
-                        if(damagelaw==1){
-                                PosPsi=max(Psi,0.);
-                                f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
-                        }
-                        else if(damagelaw==2){
-                                PosPsi=max(Psi,1.);
-                                f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
-                        }
-                        else _error_("damage law not supported");
+                        PosPsi=max(Psi,0.);
+                        f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
                 }
                 else{
                         if(equivstress==1){/* max principal stress */
@@ -311 +301 @@
                         }
                         Psi=Chi-stress_threshold;
                         NegPsi=max(-Chi,0.); /* healing only for compressive stresses */
-                        if(damagelaw==1){
-                                PosPsi=max(Psi,0.);
-                                f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
-                        }
-                        else if(damagelaw==2){
-                                PosPsi=max(Psi,1.);
-                                f[i]= c1*(pow(log10(PosPsi),c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
-                        }
-                        else _error_("damage law not supported");
+                        PosPsi=max(Psi,0.);
+                        f[i]= c1*(pow(PosPsi,c2) - healing*pow(NegPsi,c2))*pow((1./(1.-damage)),c3);
                 }
         }
         /*Add input*/
@@ -533 +516 @@
                         this->CreateDamageFInputPralong(element);
                         break;
                 case 2:
-                        this->CreateDamageFInputPralong(element);
-                        break;
-                case 3:
                         this->CreateDamageFInputExp(element);
                         break;
                 default:

../trunk-jpl/src/c/modules/Damagex/Damagex.cpp

@@ +1 @@
+/*!\file Damagex
+ * \brief: compute damage
+ */
+
+#include "./Damagex.h"
+#include "../../shared/shared.h"
+#include "../../toolkits/toolkits.h"
+
+void Damagex(FemModel* femmodel){
+
+        /*Recover Damage law Enum*/
+        int damagelaw;
+        femmodel->parameters->FindParam(&damagelaw,DamageLawEnum);
+
+        /*Calculate damage*/
+        switch(damagelaw){
+                case 0:
+                        if(VerboseModule()) _printf0_("   computing damage analytically\n");
+                        femmodel->ElementOperationx(&Element::ComputeNewDamage);
+                        break;
+                case 1:
+                case 2:
+                        if(VerboseModule()) _printf0_("   computing damage using source term in advection scheme\n");
+                        /* Damage calculated using source term in DamageEvolutionAnalysis */
+                        break;
+                default:
+                        _error_("Damage law "<<EnumToStringx(damagelaw)<<" not implemented yet");
+        }
+}

../trunk-jpl/src/c/modules/Damagex/Damagex.h

@@ +1 @@
+#ifndef _DAMAGEX_H
+#define _DAMGEX_H
+
+#include "../../classes/classes.h"
+#include "../../analyses/analyses.h"
+
+/* local prototypes: */
+void Damagex(FemModel* femmodel);
+
+#endif

../trunk-jpl/src/c/modules/modules.h

@@ -19 +19 @@
 #include "./ControlInputSetGradientx/ControlInputSetGradientx.h"
 #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h"
 #include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h"
+#include "./Damagex/Damagex.h"
 #include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
 #include "./ExpToLevelSetx/ExpToLevelSetx.h"
 #include "./ElementConnectivityx/ElementConnectivityx.h"

../trunk-jpl/src/c/Makefile.am

@@ -366 +366 @@
 endif
 if DAMAGEEVOLUTION
 issm_sources += ./analyses/DamageEvolutionAnalysis.cpp
+issm_sources += ./modules/Damagex/Damagex.cpp
 endif
 if STRESSBALANCE
 issm_sources += ./analyses/StressbalanceAnalysis.cpp

../trunk-jpl/src/c/cores/damage_core.cpp

@@ -27 +27 @@
         femmodel->parameters->FindParam(&stabilization,DamageStabilizationEnum);
 
         if(VerboseSolution()) _printf0_("   computing damage\n");
+        Damagex(femmodel); /* optionally calculate damage analytically first */
         femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
         if(stabilization==4){
                 solutionsequence_fct(femmodel);