diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 8e515fe1822b8cc2ea07a90ef36a818e89eeca05..8d76ac83bcb15523fdb6c180ecaf5c32590c40c9 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -110,6 +110,7 @@ protected:
// MatrixFunc L2_;
// MatrixFunc L3_;
+ // Elasticity tensors (in voigt notation)
FieldMatrix<double,6,6> L1_;
FieldMatrix<double,6,6> L2_;
FieldMatrix<double,6,6> L3_;
@@ -137,20 +138,20 @@ protected:
//Transformation matrix for Voigt notation
- // FieldMatrix<double,6,6> D = {{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, sqrt(2.0) , 0.0, 0.0},
- // {0.0, 0.0, 0.0, 0.0 , sqrt(2.0), 0.0},
- // {0.0, 0.0, 0.0, 0.0 , 0.0, sqrt(2.0)}
- // };
- // FieldMatrix<double,6,6> D_inv = {{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, 1.0/sqrt(2.0) , 0.0, 0.0},
- // {0.0, 0.0, 0.0, 0.0 , 1.0/sqrt(2.0), 0.0},
- // {0.0, 0.0, 0.0, 0.0 , 0.0, 1.0/sqrt(2.0)}
- // };
+ FieldMatrix<double,6,6> D_ = {{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, sqrt(2.0) , 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0 , sqrt(2.0), 0.0},
+ {0.0, 0.0, 0.0, 0.0 , 0.0, sqrt(2.0)}
+ };
+ FieldMatrix<double,6,6> D_inv = {{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, 1.0/sqrt(2.0) , 0.0, 0.0},
+ {0.0, 0.0, 0.0, 0.0 , 1.0/sqrt(2.0), 0.0},
+ {0.0, 0.0, 0.0, 0.0 , 0.0, 1.0/sqrt(2.0)}
+ };
@@ -238,12 +239,13 @@ public:
// { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[2]}
// };
-
+ // printmatrix(std::cout, D_, "D_", "--");
+ // printmatrix(std::cout, D_inv, "D_inv", "--");
// printmatrix(std::cout, L1_, "L1_", "--");
// L1_.leftmultiply(D_inv);
- // printmatrix(std::cout, L1_, "L1_", "--");
- // L1_.rightmultiply(D);
- // printmatrix(std::cout, L1_, "L1_", "--");
+ // printmatrix(std::cout, L1_, "L1_.leftmultiply(D_inv)", "--");
+ // L1_.rightmultiply(D_);
+ // printmatrix(std::cout, L1_, "L1_.rightmultiply(D_)", "--");
// Python::Module module = Python::import(materialFunctionName_);
@@ -454,16 +456,6 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
- // --- Helpers
- static FieldVector<double,6> MatrixToVoigt(const MatrixRT& G)
- {
- return {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
- }
-
- static MatrixRT VoigtToMatrix(const FieldVector<double,6>& v)
- {
- return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
- }
////////////////////////////////////////////////////////
@@ -572,26 +564,39 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
// --- Definition of transversely isotropic elasticity tensor (in voigt notation)
// - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
- // - Output: compliance Matrix
+ // - Output: inverse compliance Matrix
FieldMatrix<double,6,6> general_anisotropic(const VectorRT& framevector1,
const VectorRT& framevector2,
const VectorRT& framevector3,
const FieldMatrix<double,6,6>& C //compliance matrix
)
{
- //--- return elasticity tensor
+ //--- return elasticity tensor (in Voigt notation)
auto tmp = C;
- tmp.invert();
+ tmp.invert();
// printmatrix(std::cout, C, "C after .invert()", "--");
return tmp;
}
//-------------------------------------------------------------------------------
+ // --- Helpers
+ static FieldVector<double,6> MatrixToVoigt(const MatrixRT& G)
+ {
+ // return {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
+ return {G[0][0], G[1][1], G[2][2], 2.0*G[1][2], 2.0*G[0][2], 2.0*G[0][1]};
+
+ }
+
+ static MatrixRT VoigtToMatrix(const FieldVector<double,6>& v)
+ {
+ // return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
+ return {{v[0], v[5]/2.0, v[4]/2.0}, {v[5]/2.0, v[1], v[3]/2.0}, {v[4]/2.0, v[3]/2.0, v[2]}};
+ }
- MatrixRT applyElasticityTensor(const MatrixRT& G, const int& phase) const
+ MatrixRT applyElasticityTensor(const MatrixRT& G, const int& phase) const
{
// FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], sqrt(2.0)*G[1][2], sqrt(2.0)*G[0][2], sqrt(2.0)*G[0][1]};
// FieldVector<double,6> G_tmp = {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 902e0b882f4972fc7714e0c7f97c122438b7ffc5..8ade7d03feb08d7b4be03215d2b2638bdd928aa3 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 2 # computes all levels from first to second entry
+numLevels= 2 3 # computes all levels from first to second entry
#############################################
@@ -35,12 +35,12 @@ numLevels= 2 2 # computes all levels from first to second entry
# --- Choose material definition:
#materialFunction = "material_test"
-#materialFunction = "material_neukamm"
+materialFunction = "material_neukamm"
#materialFunction = "two_phase_material_1"
#materialFunction = "two_phase_material_2"
#materialFunction = "two_phase_material_3"
#materialFunction = "material_orthotropic"
-materialFunction = "homogeneous"
+#materialFunction = "homogeneous"
# --- Choose scale ratio gamma:
diff --git a/materials/material_neukamm.py b/materials/material_neukamm.py
index 85c128f723919481aa44d4ada2802ae44a46aeae..c654bd8b086b7123b36fb40d5b1c63c1343d0035 100644
--- a/materials/material_neukamm.py
+++ b/materials/material_neukamm.py
@@ -28,7 +28,7 @@ def indicatorFunction(x):
###############
# Wood
###############
-# def f(x):
+# def indicatorFunction(x):
# theta=0.25
# if ((abs(x[0]) < theta/2) and x[2]<0.25):
# return 1 #latewood
@@ -43,7 +43,7 @@ def indicatorFunction(x):
#########################################################################
## Notation - Parameter input :
# isotropic (Lame parameters) : [mu , lambda]
-# orthotropic : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# orthotropic : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23] # see https://en.wikipedia.org/wiki/Poisson%27s_ratio with x=1,y=2,z=3
# transversely_isotropic : [E1,E2,G12,nu12,nu23]
# general_anisotropic : full compliance matrix C
######################################################################
diff --git a/materials/material_orthotropic.py b/materials/material_orthotropic.py
index 6e68cb6451034044f23c66b8c132a96d138fdf4d..ffaf3901ff1af4ce75a564180c7b064efa4897b2 100644
--- a/materials/material_orthotropic.py
+++ b/materials/material_orthotropic.py
@@ -24,7 +24,7 @@ def indicatorFunction(x):
#########################################################################
## Notation - Parameter input :
# isotropic (Lame parameters) : [mu , lambda]
-# orthotropic : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# orthotropic : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23] # see https://en.wikipedia.org/wiki/Poisson%27s_ratio with x=1,y=2,z=3
# transversely_isotropic : [E1,E2,G12,nu12,nu23]
# general_anisotropic : full compliance matrix C
######################################################################
diff --git a/materials/material_test.py b/materials/material_test.py
index 7d59bae264b5b49604cb6165fe30fd566fbeea07..b3f0737dcdef5f81fd7383f5833bfa2f43494acf 100644
--- a/materials/material_test.py
+++ b/materials/material_test.py
@@ -16,7 +16,7 @@ def indicatorFunction(x):
elif (x[1]<-0.5+theta and x[2]>0.5-theta):
return 2 #Phase2
else :
- return 0 #Phase3
+ return 3 #Phase3
########### Options for material phases: #################################
@@ -24,7 +24,7 @@ def indicatorFunction(x):
#########################################################################
## Notation - Parameter input :
# isotropic (Lame parameters) : [mu , lambda]
-# orthotropic : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23]
+# orthotropic : [E1,E2,E3,G12,G23,G31,nu12,nu13,nu23] # see https://en.wikipedia.org/wiki/Poisson%27s_ratio with x=1,y=2,z=3
# transversely_isotropic : [E1,E2,G12,nu12,nu23]
# general_anisotropic : full compliance matrix C
######################################################################
@@ -36,32 +36,40 @@ Phases=3
#--- Define different material phases:
#- PHASE 1
-phase1_type="isotropic"
-materialParameters_phase1 = [80, 80]
+# phase1_type="isotropic"
+# materialParameters_phase1 = [80, 80]
# phase1_type="orthotropic"
# materialParameters_phase1 = [11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37] # walnut parameters (values for compliance matrix)
# # materialParameters_phase1 = [10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31] # Norway spruce parameters (values for compliance matrix)
+phase1_type="general_anisotropic"
+materialParameters_phase1 = np.array([[1.0, 8.0, 16.0, 16.0 , 8.0, 8.0],
+ [8.0, 1.0, 16.0, 16.0 , 8.0, 8.0],
+ [8.0, 8.0, 1.0, 16.0 , 8.0, 8.0],
+ [8.0, 8.0, 16.0, math.sqrt(2), 8.0, 8.0],
+ [8.0, 8.0, 16.0, 16.0 , 8.0, 8.0],
+ [8.0, 8.0, 16.0, 16.0 , 8.0, 1.0]])
+
#- PHASE 2
-# phase2_type="transversely_isotropic"
-# materialParameters_phase2 = [11.2e3,1190,960,0.63 ,0.37]
+phase2_type="transversely_isotropic"
+materialParameters_phase2 = [11.2e3,1190,960,0.63 ,0.37]
-phase2_type="isotropic"
-materialParameters_phase2 = [80, 80]
+# phase2_type="isotropic"
+# materialParameters_phase2 = [80, 80]
#- PHASE 3
-phase3_type="isotropic"
-materialParameters_phase3 = [60, 25]
+# 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]])
+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