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

Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp     (revision 19387)
+++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp     (revision 19388)
@@ -129,7 +129,7 @@
        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,7 +139,6 @@
        element->FindParam(&epsf,DamageC1Enum);
        element->FindParam(&stress_threshold,DamageStressThresholdEnum);
        element->FindParam(&domaintype,DomainTypeEnum);
-       element->FindParam(&damagelaw,DamageLawEnum);
 
        /*Compute stress tensor: */
        element->ComputeStrainRate();
@@ -201,7 +200,7 @@
        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,7 +213,6 @@
        element->FindParam(&healing,DamageHealingEnum);
        element->FindParam(&stress_threshold,DamageStressThresholdEnum);
        element->FindParam(&domaintype,DomainTypeEnum);
-       element->FindParam(&damagelaw,DamageLawEnum);
 
        /*Get problem dimension*/
        switch(domaintype){
@@ -290,16 +288,8 @@
                        }
                        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,15 +301,8 @@
                        }
                        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,9 +516,6 @@
                        this->CreateDamageFInputPralong(element);
                        break;
                case 2:
-                       this->CreateDamageFInputPralong(element);
-                       break;
-               case 3:
                        this->CreateDamageFInputExp(element);
                        break;
                default:
Index: ../trunk-jpl/src/c/modules/Damagex/Damagex.cpp
===================================================================
--- ../trunk-jpl/src/c/modules/Damagex/Damagex.cpp      (revision 0)
+++ ../trunk-jpl/src/c/modules/Damagex/Damagex.cpp      (revision 19388)
@@ -0,0 +1,29 @@
+/*!\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");
+       }
+}
Index: ../trunk-jpl/src/c/modules/Damagex/Damagex.h
===================================================================
--- ../trunk-jpl/src/c/modules/Damagex/Damagex.h        (revision 0)
+++ ../trunk-jpl/src/c/modules/Damagex/Damagex.h        (revision 19388)
@@ -0,0 +1,10 @@
+#ifndef _DAMAGEX_H
+#define _DAMGEX_H
+
+#include "../../classes/classes.h"
+#include "../../analyses/analyses.h"
+
+/* local prototypes: */
+void Damagex(FemModel* femmodel);
+
+#endif

Property changes on: ../trunk-jpl/src/c/modules/Damagex
___________________________________________________________________
Added: svn:ignore
   + *.deps
*.dirstamp


Index: ../trunk-jpl/src/c/modules/modules.h
===================================================================
--- ../trunk-jpl/src/c/modules/modules.h        (revision 19387)
+++ ../trunk-jpl/src/c/modules/modules.h        (revision 19388)
@@ -19,6 +19,7 @@
 #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"
Index: ../trunk-jpl/src/c/Makefile.am
===================================================================
--- ../trunk-jpl/src/c/Makefile.am      (revision 19387)
+++ ../trunk-jpl/src/c/Makefile.am      (revision 19388)
@@ -366,6 +366,7 @@
 endif
 if DAMAGEEVOLUTION
 issm_sources += ./analyses/DamageEvolutionAnalysis.cpp
+issm_sources += ./modules/Damagex/Damagex.cpp
 endif
 if STRESSBALANCE
 issm_sources += ./analyses/StressbalanceAnalysis.cpp
Index: ../trunk-jpl/src/c/cores/damage_core.cpp
===================================================================
--- ../trunk-jpl/src/c/cores/damage_core.cpp    (revision 19387)
+++ ../trunk-jpl/src/c/cores/damage_core.cpp    (revision 19388)
@@ -27,6 +27,7 @@
        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);