Update + Add more materials

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;          
             }

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."); 
+     
   }
 
 

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
 
 
 

materials/homogeneous.py

+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]]

materials/material_neukamm.py

 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]]]

materials/material_neukamm_old.py

+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]]

materials/material_orthotropic.py

+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]]]

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):

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]]

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]]

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]]

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