From cf7da5e4dcf1dbeabc06913016aa51cd3161600e Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 31 Aug 2022 20:37:04 +0200
Subject: [PATCH] Update + Add more materials

---
 .../EffectiveQuantitiesComputer.hh            |   8 +-
 dune/microstructure/prestrainedMaterial.hh    |   8 +-
 inputs/cellsolver.parset                      |  16 ++--
 .../material_neukamm.cpython-38.pyc           | Bin 821 -> 1179 bytes
 .../__pycache__/material_test.cpython-38.pyc  | Bin 1177 -> 1176 bytes
 .../two_phase_material_1.cpython-38.pyc       | Bin 324 -> 788 bytes
 .../two_phase_material_2.cpython-38.pyc       | Bin 331 -> 783 bytes
 materials/homogeneous.py                      |  35 +++++++++
 materials/material_neukamm.py                 |  73 ++++++++++++------
 materials/material_neukamm_old.py             |  50 ++++++++++++
 materials/material_orthotropic.py             |  71 +++++++++++++++++
 materials/material_test.py                    |   2 +-
 materials/two_phase_material_1.py             |  35 ++++++++-
 materials/two_phase_material_2.py             |  35 ++++++++-
 materials/two_phase_material_3.py             |  31 ++++++++
 src/Cell-Problem.cc                           |   1 +
 .../material_backup.py                        |   0
 17 files changed, 324 insertions(+), 41 deletions(-)
 create mode 100644 materials/homogeneous.py
 create mode 100644 materials/material_neukamm_old.py
 create mode 100644 materials/material_orthotropic.py
 rename {materials => src/deprecated_code/31-8-22_prestrainedMaterials}/material_backup.py (100%)

diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index 0d56aade..9329f163 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -60,7 +60,7 @@ public:
     { 
 
         // prestrain_ = material_.getPrestrainFunction();
-        std::cout << "starting effective quantities computation..." << std::endl;
+        
 
     } 
 
@@ -79,6 +79,8 @@ public:
     void computeEffectiveQuantities()
     {
 
+        std::cout << "starting effective quantities computation..." << std::endl;
+
         // Get everything.. better TODO: with Inheritance?
         // auto test = correctorComputer_.getLoad_alpha1();
         // auto phiContainer = correctorComputer_.getPhicontainer();
@@ -179,6 +181,7 @@ public:
                     elementEnergy += energyDensity * quadPoint.weight() * integrationElement;      // quad[quadPoint].weight() ???
                     if (b==0)
                     {
+                        // std::cout << "WORKS TILL HERE" << std::endl;
                         // auto prestrainPhaseValue_old = prestrain_(localIndicatorFunction(quadPos));
                         auto prestrainPhaseValue = material_.prestrain(localIndicatorFunction(quadPos),element.geometry().global(quadPos));
                         // if (localIndicatorFunction(quadPos) != 3)
@@ -291,7 +294,8 @@ public:
                 const auto& quadPos = quadPoint.position();
                 const double integrationElement = geometry.integrationElement(quadPos);
 
-                double energyDensity= scalarProduct(material_.ElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
+                double energyDensity= scalarProduct(material_.applyElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
+                // double energyDensity= scalarProduct(material_.ElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
                 // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), matrixFieldA(quadPos), matrixFieldB(quadPos));
                 elementEnergy += energyDensity * quadPoint.weight() * integrationElement;          
             }
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 4a020ed9..8e515fe1 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -30,9 +30,6 @@ using std::make_shared;
 
 
 
-
-
-
 //   // ----------------------------------------------------------------------------------
 //   template<class Domain>
 //   int indicatorFunction_material(const Domain& x)
@@ -635,8 +632,11 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
       return prestrain1_(x);
     else if (phase == 2)
       return prestrain2_(x);
-    else 
+    else if (phase ==3)
       return prestrain3_(x);
+    else 
+      DUNE_THROW(Exception, "Phase number not implemented."); 
+     
   }
 
 
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 9ce849f0..902e0b88 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -26,7 +26,7 @@ outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
 ## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
 #----------------------------------------------------
 
-numLevels= 2 3      # computes all levels from first to second entry
+numLevels= 2 2      # computes all levels from first to second entry
 
 
 #############################################
@@ -34,12 +34,14 @@ numLevels= 2 3      # computes all levels from first to second entry
 #############################################
 
 # --- Choose material definition:
+#materialFunction  = "material_test"
+#materialFunction  = "material_neukamm"
 #materialFunction = "two_phase_material_1"
 #materialFunction = "two_phase_material_2"
-#materialFunction  = "material_neukamm"
-materialFunction  = "material"
-#materialFunction = "material_2"
-#materialFunction = "homogeneous"
+#materialFunction = "two_phase_material_3"
+#materialFunction = "material_orthotropic"
+materialFunction = "homogeneous"
+
 
 # --- Choose scale ratio gamma:
 gamma=1.0
@@ -47,8 +49,8 @@ gamma=1.0
 
 
 
-#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
-geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
+#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/materials
+geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/materials
 
 
 
diff --git a/materials/__pycache__/material_neukamm.cpython-38.pyc b/materials/__pycache__/material_neukamm.cpython-38.pyc
index f077d3db1ae6d8c15bcd7f2a1addeacdd7531020..9d7421aeea4ce841358b2a6f70980f04d0e1084f 100644
GIT binary patch
literal 1179
zcmWIL<>g{vU|_hqm_M1BgMr~Oh=Yuo7#J8F7#J9eOBfg!QW#Pga~PsPG*b>^2#99P
zVTxh`i!(<tgK3s1mK5d`mKKI6))dw#wiLD~_7wIM_H?Ewjuf^iPAHozg)NFZg(HPC
zohgbZg)4=-g&~SJg(ro#g&~SBg)fD_g&~SRg)x{xQ{W}YJ${;uw-~jO89{6)W@cbu
z0AXj4^CCdbW2|LNVX9@SVMt+?WJqDDVH9V`W+>9BVMt*`=BF?PGt@AoFd>NrK-I9-
zFiJ3FvlIoSum&@vFb6ZFFa$Gbviqe&?7U$A66As+5CIi>u>U0}j9@(bl?<B9w-_s;
zSW7ZeOA@2l(h`$P@{5XC85kJA#IGRzjQreG{p_5?(qes=)Z*-t`~rQK(!5mNkksOm
z%)E5{6c8sjGr1_gxTL5wxumoxRX;beB(*3rF{c>LiqA_e%}&hC)hnpHC779)l9`-X
zl3(OjnwMOXnV%;F3KTw&s~K4sSr`={n5l?^fq~&ADAb}jGmG;}it-CGlV1ja1Z_Zs
zBm)BjN+5%L3UWI*kY&JuoWhjCl+K*Ol)?z5nNk?jnbVnr8CEj;6@knKhY#4dD;bN}
zU_J*66tOcfFx(O-C`v6ZDN4-Di!aDXEKW6a0XY{Gjtq=N{1myNh!<`j$omjWK|a84
zvynTJ%>u;O%n0&r3L^-nGlNJFALM&YCO?G7i-@t@*bB*WK?Vi}O@Ui1xrrqiQ4$4}
zB^mj7@wtg5MVS@x`30#(iJ)jLj$%tLsVqn>zQvSZe2clbvKSnrQEUOAU@VT}h6HCk
zNKKR^EYSoc7A5AUmZTOHLxUTt&WM0IBd9uK0_u#xaaIHhz9Nt}kOCbM@*Fm~`6;D2
gsdk_wRV>HAz`(%>g}fY!EX*v795Bqv!OX!70K+2t#sB~S

delta 378
zcmbQuxs{DCl$V!_fq{Wx`Y+C8b(V>IG8Q@v3=Am@DU3M`Q4A?eDa<VlQH&`pDXc9F
zQA{apDeNr_QOqeEDV!|~Q7kEp!3>&QFF^+QO<cc=n}vaap$J6$3Z7iTXv&$MnxC6m
zQk0oGc_X8kG-H|q0|P@bA4oYP3nL4o0t7P^u}>Ccax`TE>1G5GAR|^X7O}$EU^9!@
z7#J9CF(nylflL7@W?(Ghnk>j<t6#(cSIG%yfeeA^G15iS!#!D$NoBGzvr;)&33eUE
z21q)17#J8dS#PoACYEFrv4d;?1zHg|hy{{^xDRX^m;k$!!zMRBr8Fni4&<d`kPA8F
KIk-5OIhX-e@H;2~

diff --git a/materials/__pycache__/material_test.cpython-38.pyc b/materials/__pycache__/material_test.cpython-38.pyc
index 1dfcaeea6e9e18505d515bf97694734c7d04784c..5185825128d19bce06d794b7260599da858ed7ed 100644
GIT binary patch
delta 52
zcmbQqIfIiol$V!_fq{XcV-bJyoQ=F`jEw%1OBhW#auZ8Zi!u{)Hg96Q$Hd9Z!pOnL
I!OX!709`f=jsO4v

delta 53
zcmbQiIg^t&l$V!_fq{X+dl7$f(MH}hM#g~2C5)z=>8bg-sU=03shc-4-eclqVPWLp
J<zVJu1^{0k44MD{

diff --git a/materials/__pycache__/two_phase_material_1.cpython-38.pyc b/materials/__pycache__/two_phase_material_1.cpython-38.pyc
index c9763e8459ecb9c81f61150564e7b902e4558747..78b1dd2bd3091a8f4709ecd09742ae4c3d6ac0a7 100644
GIT binary patch
literal 788
zcmWIL<>g{vU|@K(oIlxum4V?gh=Yuo7#J8F7#J9eK|~5e3S$mK6hjJA3UdoX6k`fY
z6jKUo6mtq&3R^l;6iW(g6f2a?mcpLG(ZUeLp2C^J)xr?Pk-`|vpvnCbWVT;2BZz`x
z5Sy8Sfx#JMiwpw;LkU9-LkgoLLo;IvQ#MPHND6Z>LkdeULkdGMgC?tA5lBNiMDqpv
zmmuCtkSjGAZ!sq(6|ZEx#aOYDp@^M<f#Fw}enx(7s(yA(Vrj9yOKNd;Nq&L8OKD!J
zZb)izNoHQUehP?_o0(jcUtCgDnp{#^l&YVbSdv<lnV3_oUs9eQUyzYloEi@ki#OCO
zsJtbZnU|88oLG`y<W`!OT#}ie#{qIZ8^~2m0*ow7MVt%_3`HPsMR8^p=a&@a7i1>C
z3}9ej!0<aOvfon}Q<zei(wS43KrArLl){+KoX#A~u#(xY2;@<40D%2h1kwv(fxKJ9
z!N9<9OQ4`AwYa1xF*6V1bVCi0(?A*-7>jr)azPO{$X8&ih_Ts78_8x~1_lOA&RZ<G
zi6t4}(1~IT0EKLE6gN184dY8H3sR#bVF4YGSd^HXT9R5+4D|+7oe=?bMqob{fviIc
o8L*qd!NXybo1apelWGS_IK=`C3=ACX96~J2EQ}l=%*MeC0Ked-i~s-t

delta 214
zcmbQjc7%y9l$V!_fq{YHRmHI+CB}(-GEyQ83=Am@DU3M`Q4A?eDa<VlQH&{!!3>%#
z6ALVz8G{*87(g7qB1Q%VhIEi#28IjvFF`y_##_vZNyVCsw-_r{G8C~gFfjb`(9g)v
zP1P?hNli=4({~B=bDhK}%9);;pPO1zl$knNmQl)wF^z+PfuWcUWEhhGBMVazI|BoQ
nCetmJ+{BU$Fb8ZH5&=?si^C>2KczG$)edBNF%ts=0|z?*Pp>QH

diff --git a/materials/__pycache__/two_phase_material_2.cpython-38.pyc b/materials/__pycache__/two_phase_material_2.cpython-38.pyc
index de9026c11a1592fe0b109938ac6617277d4588af..86b1c59d67e8cc248fd7ae49f0d16bdafe0fc4fc 100644
GIT binary patch
literal 783
zcmWIL<>g{vU|{H4!Jllw%E0g##6iYP3=9ko3=9m#AR>h!g)xUAiXnw5g}H?xiZO*H
ziYbLPiaCWXg)N;aiY0|LiWSOcOJPspXkmzAPvK1AYGH`tNMQ_S(ByszGTSei5kx^T
zh|SEvz~BtBMTUWap@gA^A%#(rp_ws-DVwE81f(;CIhY}ZA(%mv#qT9ZLpntB1^bsE
z2Wm3jVoppdUdeciv0^1d5jz6|!>=&?jQreG{p_5?(qes=)Z*-t`~rQK(!5mNkksOm
z%)E5{6c8sjGr1_gxTL5wxumoxRX;beB(*3rF{fC+q&z>qAS1CjH6A7wZ=_dHc}p-e
zFC{ZMu_V99tu!yWBr`vc1LSZvkb{^67+IK#I2jliia<V!;>;}0FDc3|$V`42z`($O
z;cZrAZ>KP(Fr_f1Gp8_tSYVndg)yBuojI6cC9_`<$fHFdzk~f(1kwv(fxKJ9!N9<9
zOQ4`AwYa1xF*6V1bVCi0(?A*-7>jr)azPO{$X8&ih_Ts78_8x~1_lOA&RZ<Gi6t4}
z(1~IT0EKLE6gN184dY8H3sR#bVF4YGSd^HXT9R5+4D|+7oe=?bMqob{fviIc8L*qd
k!NXybo1apelWGSF#bN;l1_lmx4j~q17Df&bX5(N603B|nzW@LL

delta 236
zcmeBYJI%xw%FD~ez`(#zl5-?UiE$#IjFboi149Z!3S$mK6hjJA3UdoX6k`fwFoPz`
z#Dc*3B1Q%V22I9W%!x_GnvAy?D^@ZTu`)0){PNe&$j?pH&(29KE!KBQEzT~<FVJ@h
z^>fuvDa}jO&CN_M$}cV{DorjaElSl-PtDIwEh)-OE!Hn7&yO$2NGwi`&rK{zEy_&H
zi8s<KsJz9P#=*eAP|OB$9g_ef3sVt01A{q(CetmJ+{BU$5D#KC*nS8BQhbZUCO1E&
RG$+*#<kVs&1_lNWb^!XsJI(+A

diff --git a/materials/homogeneous.py b/materials/homogeneous.py
new file mode 100644
index 00000000..7369270d
--- /dev/null
+++ b/materials/homogeneous.py
@@ -0,0 +1,35 @@
+import math
+
+
+#Indicator function that determines both phases
+# x[0] : x-component
+# x[1] : y-component
+# x[2] : z-component
+def indicatorFunction(x):
+    # --- replace with your definition of indicatorFunction:
+    return 1    #Phase1
+
+
+
+########### Options for material phases: #################################
+#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+#########################################################################
+## Notation - Parameter input :
+# isotropic (Lame parameters) : [mu , lambda]
+#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# general_anisotropic         : full compliance matrix C
+######################################################################
+# --- Number of material phases
+Phases=1
+#--- Define different material phases:
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [80, 80]
+
+
+
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
diff --git a/materials/material_neukamm.py b/materials/material_neukamm.py
index 42440c6d..85c128f7 100644
--- a/materials/material_neukamm.py
+++ b/materials/material_neukamm.py
@@ -1,11 +1,30 @@
 import math
+from python_matrix_operations import *
+import ctypes
+import os
+import sys
+#---------------------------------------------------------------
 
 
-#Indicator function that determines both phases
+
+#--- define indicator function for material phases
 # x[0] : y1-component
 # x[1] : y2-component
 # x[2] : x3-component
-#    --- replace with your definition of indicatorFunction:
+#To indicate phases return either : 1 / 2 / 3
+###############
+# Cross
+###############
+def indicatorFunction(x):
+    theta=0.25
+    factor=1
+    if (x[0] <-0.5+theta and x[2]<-0.5+theta):
+        return 1    #Phase1
+    elif (x[1]<-0.5+theta and x[2]>0.5-theta):
+        return 2    #Phase2
+    else :
+        return 3   #Phase3
+
 ###############
 # Wood
 ###############
@@ -18,32 +37,40 @@ import math
     # else :
     #     return 0    #Phase3
 
-# def b1(x):
-#     return [[.5, 0, 0], [0,1,0], [0,0,0]]
 
-# def b2(x):
-#     return [[.4, 0, 0], [0,.4,0], [0,0,0]]
+########### Options for material phases: #################################
+#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+#########################################################################
+## Notation - Parameter input :
+# isotropic (Lame parameters) : [mu , lambda]
+#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# general_anisotropic         : full compliance matrix C
+######################################################################
+# --- Number of material phases
+Phases=3
+#--- Define different material phases:
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [80, 80]
+
+#- PHASE 2
+phase2_type="isotropic"
+materialParameters_phase2 = [80, 80]
+
+#- PHASE 3
+phase3_type="isotropic"
+materialParameters_phase3 = [60, 25]
 
-# def b3(x):
-#     return [[0, 0, 0], [0,0,0], [0,0,0]]
-###############
-# Cross
-###############
-def f(x):
-    theta=0.25
-    factor=1
-    if (x[0] <-1/2+theta and x[2]<-1/2+theta):
-        return 1    #Phase1
-    elif (x[1]< -1/2+theta and x[2]>1/2-theta):
-        return 2    #Phase2
-    else :
-        return 0    #Phase3
 
-def b1(x):
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
     return [[1, 0, 0], [0,1,0], [0,0,1]]
 
-def b2(x):
+def prestrain_phase2(x):
     return [[1, 0, 0], [0,1,0], [0,0,1]]
 
-def b3(x):
+def prestrain_phase3(x):
     return [[0, 0, 0], [0,0,0], [0,0,0]]
+    # return [[x[0],0 ,0 ], [0,0,x[1]], [0,0,x[2]]]
diff --git a/materials/material_neukamm_old.py b/materials/material_neukamm_old.py
new file mode 100644
index 00000000..790e3a43
--- /dev/null
+++ b/materials/material_neukamm_old.py
@@ -0,0 +1,50 @@
+import math
+
+# DEPRECATED!!!! just for reference
+
+#Indicator function that determines both phases
+# x[0] : y1-component
+# x[1] : y2-component
+# x[2] : x3-component
+#    --- replace with your definition of indicatorFunction:
+###############
+# Wood
+###############
+# def f(x):
+#     theta=0.25
+    # if ((abs(x[0]) < theta/2) and x[2]<0.25):
+    #     return 1    #latewood
+    # elif ((abs(x[0]) > theta/2) and x[2]<0.25):
+    #     return 2    #earlywood
+    # else :
+    #     return 0    #Phase3
+
+# def b1(x):
+#     return [[.5, 0, 0], [0,1,0], [0,0,0]]
+
+# def b2(x):
+#     return [[.4, 0, 0], [0,.4,0], [0,0,0]]
+
+# def b3(x):
+#     return [[0, 0, 0], [0,0,0], [0,0,0]]
+###############
+# Cross
+###############
+def f(x):
+    theta=0.25
+    factor=1
+    if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+        return 1    #Phase1
+    elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+        return 2    #Phase2
+    else :
+        return 0    #Phase3
+
+def b1(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def b2(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def b3(x):
+    return [[0, 0, 0], [0,0,0], [0,0,0]]
diff --git a/materials/material_orthotropic.py b/materials/material_orthotropic.py
new file mode 100644
index 00000000..6e68cb64
--- /dev/null
+++ b/materials/material_orthotropic.py
@@ -0,0 +1,71 @@
+import math
+from python_matrix_operations import *
+import ctypes
+import os
+import sys
+#---------------------------------------------------------------
+
+
+
+#--- define indicator function
+def indicatorFunction(x):
+    theta=0.25
+    factor=1
+    if (x[0] <-0.5+theta and x[2]<-0.5+theta):
+        return 1    #Phase1
+    elif (x[1]<-0.5+theta and x[2]>0.5-theta):
+        return 2    #Phase2
+    else :
+        return 3    #Phase3
+
+
+########### Options for material phases: #################################
+#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+#########################################################################
+## Notation - Parameter input :
+# isotropic (Lame parameters) : [mu , lambda]
+#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# general_anisotropic         : full compliance matrix C
+######################################################################
+
+# --- Number of material phases
+Phases=3
+
+
+#--- Define different material phases:
+
+#- PHASE 1
+phase1_type="orthotropic"
+materialParameters_phase1 = [11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37]    # walnut parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+
+#- PHASE 2
+phase2_type="orthotropic"
+materialParameters_phase2 = [10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31]   # Norway spruce parameters (values for compliance matrix) see [Dimwoodie; Timber its nature and behavior p.109]
+
+#- PHASE 3
+phase3_type="isotropic"
+materialParameters_phase3 = [60, 25]
+
+
+#--- for general anisotopic material the compliance matrix is required:
+# phase3_type="general_anisotropic"
+# materialParameters_phase3 = np.array([[1.0,     0.0,     0.0,   0.0 ,         0.0,   0.0],
+#                                       [0.0,     1.0,     0.0,   0.0 ,         0.0,   0.0],
+#                                       [0.0,     0.0,     1.0,   0.0 ,         0.0,   0.0],
+#                                       [0.0,     0.0,     0.0,   math.sqrt(2), 0.0,   0.0],
+#                                       [0.0,     0.0,     0.0,   0.0 ,         1.0,   0.0],
+#                                       [0.0,     0.0,     0.0,   0.0 ,         0.0,   1.0]])
+
+
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def prestrain_phase2(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def prestrain_phase3(x):
+    return [[0, 0, 0], [0,0,0], [0,0,0]]
+    # return [[x[0],0 ,0 ], [0,0,x[1]], [0,0,x[2]]]
diff --git a/materials/material_test.py b/materials/material_test.py
index fed02ba5..7d59bae2 100644
--- a/materials/material_test.py
+++ b/materials/material_test.py
@@ -5,7 +5,7 @@ import os
 import sys
 #---------------------------------------------------------------
 
-
+#To indicate phases return either : 1 / 2 / 3
 
 #--- define indicator function
 def indicatorFunction(x):
diff --git a/materials/two_phase_material_1.py b/materials/two_phase_material_1.py
index 8c3418f4..b09bf3c9 100644
--- a/materials/two_phase_material_1.py
+++ b/materials/two_phase_material_1.py
@@ -5,9 +5,40 @@ import math
 # x[0] : x-component
 # x[1] : y-component
 # x[2] : z-component
-def f(x):
+def indicatorFunction(x):
     # --- replace with your definition of indicatorFunction:
     if (abs(x[0]) > 0.25):
         return 1    #Phase1
     else :
-        return 0    #Phase2
+        return 2    #Phase2
+
+
+########### Options for material phases: #################################
+#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+#########################################################################
+## Notation - Parameter input :
+# isotropic (Lame parameters) : [mu , lambda]
+#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# general_anisotropic         : full compliance matrix C
+######################################################################
+# --- Number of material phases
+Phases=2
+#--- Define different material phases:
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [80, 80]
+
+#- PHASE 2
+phase2_type="isotropic"
+materialParameters_phase2 = [80, 80]
+
+
+
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def prestrain_phase2(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
diff --git a/materials/two_phase_material_2.py b/materials/two_phase_material_2.py
index b643b490..b26dd94b 100644
--- a/materials/two_phase_material_2.py
+++ b/materials/two_phase_material_2.py
@@ -5,9 +5,40 @@ import math
 # x[0] : x-component
 # x[1] : y-component
 # x[2] : z-component
-def f(x):
+def indicatorFunction(x):
     # --- replace with your definition of indicatorFunction:
     if (abs(x[1]) > 0.25):
         return 1    #Phase1
     else :
-        return 0    #Phase2
+        return 2    #Phase2
+
+
+########### Options for material phases: #################################
+#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+#########################################################################
+## Notation - Parameter input :
+# isotropic (Lame parameters) : [mu , lambda]
+#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# general_anisotropic         : full compliance matrix C
+######################################################################
+# --- Number of material phases
+Phases=2
+#--- Define different material phases:
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [80, 80]
+
+#- PHASE 2
+phase2_type="isotropic"
+materialParameters_phase2 = [80, 80]
+
+
+
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def prestrain_phase2(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
diff --git a/materials/two_phase_material_3.py b/materials/two_phase_material_3.py
index fea31cb1..f087f5d4 100644
--- a/materials/two_phase_material_3.py
+++ b/materials/two_phase_material_3.py
@@ -11,3 +11,34 @@ def f(x):
         return 1    #Phase1
     else :
         return 0    #Phase2
+
+
+########### Options for material phases: #################################
+#     1. "isotropic"     2. "orthotropic"      3. "transversely_isotropic"   4. "general_anisotropic"
+#########################################################################
+## Notation - Parameter input :
+# isotropic (Lame parameters) : [mu , lambda]
+#         orthotropic         : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# transversely_isotropic      : [E1,E2,G12,nu12,nu23]
+# general_anisotropic         : full compliance matrix C
+######################################################################
+# --- Number of material phases
+Phases=2
+#--- Define different material phases:
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [80, 80]
+
+#- PHASE 2
+phase2_type="isotropic"
+materialParameters_phase2 = [80, 80]
+
+
+
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def prestrain_phase2(x):
+    return [[1, 0, 0], [0,1,0], [0,0,1]]
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index af299a97..9899928e 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -235,6 +235,7 @@ int main(int argc, char *argv[])
 
     // --- get scale ratio 
     double gamma = parameterSet.get<double>("gamma",1.0); 
+    std::cout << "scale ratio (gamma) set to : " << gamma << std::endl;
 
     //------------------------------------------------------------------------------------------------
     //--- Compute Correctors
diff --git a/materials/material_backup.py b/src/deprecated_code/31-8-22_prestrainedMaterials/material_backup.py
similarity index 100%
rename from materials/material_backup.py
rename to src/deprecated_code/31-8-22_prestrainedMaterials/material_backup.py
-- 
GitLab