From 256231c24a38730d307da19fdbc0a751a7a6e482 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Fri, 29 Sep 2023 15:38:26 +0200
Subject: [PATCH] Increase feasible number of material phases to 4 and add
parametrized_laminate
---
dune/microstructure/prestrainedMaterial.hh | 198 +++++++++++++++------
inputs/cellsolver.parset | 11 +-
materials/parametrized_laminate.py | 68 +++++++
3 files changed, 216 insertions(+), 61 deletions(-)
create mode 100644 materials/parametrized_laminate.py
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 904b1e58..3634e4fe 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -113,12 +113,12 @@ protected:
Dune::FieldMatrix<double,6,6> L1_;
Dune::FieldMatrix<double,6,6> L2_;
Dune::FieldMatrix<double,6,6> L3_;
- Dune::FieldMatrix<double,6,6> L4_; //not used yet
+ Dune::FieldMatrix<double,6,6> L4_;
Func2Tensor prestrain1_;
Func2Tensor prestrain2_;
Func2Tensor prestrain3_;
- Func2Tensor prestrain4_; //not used yet
+ Func2Tensor prestrain4_;
Python::Module module_;
@@ -138,19 +138,19 @@ protected:
//Transformation matrix for Voigt notation
Dune::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)}
- };
+ {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)}
+ };
Dune::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)}
- };
+ {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)}
+ };
@@ -287,6 +287,11 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
materialParameters_string= "materialParameters_phase3";
break;
}
+ case 4:
+ {
+ materialParameters_string= "materialParameters_phase4";
+ break;
+ }
default:
DUNE_THROW(Dune::Exception, "Phase number not implemented.");
break;
@@ -296,6 +301,7 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
{
Dune::FieldVector<double,2> materialParameters(0);
module.get(materialParameters_string).toC<Dune::FieldVector<double,2>>(materialParameters);
+ std::cout << "Lame-Parameters (mu,lambda): (" << materialParameters[0] << "," << materialParameters[1] << ") " << std::endl;
return isotropic(materialParameters[0],materialParameters[1]);
}
if(phaseType == "orthotropic") //TODO SetupPHASE (phase_type)
@@ -394,24 +400,29 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
prestrain3_ = Python::make_function<MatrixRT>(module.get("prestrain_phase3"));
- // if (phase1_type == "isotropic" ) //TODO SetupPHASE (phase_type)
- // {
- // FieldVector<double,2> materialParameters_phase1(0);
- // module.get("materialParameters_phase1").toC<FieldVector<double,2>>(materialParameters_phase1);
- // L1_ = isotropic(materialParameters_phase1[0],materialParameters_phase1[1]);
- // }
- // if (phase2_type == "isotropic" ) //TODO SetupPHASE (phase_type)
- // {
- // FieldVector<double,2> materialParameters_phase2(0);
- // module.get("materialParameters_phase2").toC<FieldVector<double,2>>(materialParameters_phase2);
- // L2_ = isotropic(materialParameters_phase2[0],materialParameters_phase2[1]);
- // }
- // if (phase3_type == "isotropic" ) //TODO SetupPHASE (phase_type)
- // {
- // FieldVector<double,2> materialParameters_phase3(0);
- // module.get("materialParameters_phase3").toC<FieldVector<double,2>>(materialParameters_phase3);
- // L3_ = isotropic(materialParameters_phase3[0],materialParameters_phase3[1]);
- // }
+ break;
+ }
+ case 4:
+ {
+ std::string phase1_type;
+ module.get("phase1_type").toC<std::string>(phase1_type);
+ std::string phase2_type;
+ module.get("phase2_type").toC<std::string>(phase2_type);
+ std::string phase3_type;
+ module.get("phase3_type").toC<std::string>(phase3_type);
+ std::string phase4_type;
+ module.get("phase4_type").toC<std::string>(phase4_type);
+
+ L1_ = setupPhase(phase1_type, module,1);
+ L2_ = setupPhase(phase2_type, module,2);
+ L3_ = setupPhase(phase3_type, module,3);
+ L4_ = setupPhase(phase4_type, module,4);
+
+ prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
+ prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
+ prestrain3_ = Python::make_function<MatrixRT>(module.get("prestrain_phase3"));
+ prestrain4_ = Python::make_function<MatrixRT>(module.get("prestrain_phase4"));
+
break;
}
default:
@@ -507,7 +518,7 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
{ 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 1.0/G_12}
};
-
+ printmatrix(std::cout, C, "(Orthotropic) Compliance matrix used:", "--");
// printmatrix(std::cout, C, "C", "--");
//--- return elasticity tensor
@@ -593,6 +604,32 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
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]}};
}
+ static Dune::FieldVector<double,6> strainToVoigt(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 Dune::FieldVector<double,6> stressToVoigt(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], G[1][2], G[0][2], G[0][1]};
+
+ }
+
+ static MatrixRT VoigtToStress(const Dune::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], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
+ }
+
+ static MatrixRT VoigtToStrain(const Dune::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
@@ -603,25 +640,57 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
Dune::FieldVector<double,6> out(0);
// printvector(std::cout, G_tmp, "G_tmp", "--");
- if (phase == 1)
- {
- // printmatrix(std::cout, L1_, "L1_", "--");
- L1_.mv(G_tmp,out);
- }
- else if (phase == 2)
+ // if (phase == 1)
+ // {
+ // // printmatrix(std::cout, L1_, "L1_", "--");
+ // L1_.mv(G_tmp,out);
+ // }
+ // else if (phase == 2)
+ // {
+ // // printmatrix(std::cout, L2_, "L2_", "--");
+ // L2_.mv(G_tmp,out);
+ // }
+ // else
+ // {
+ // // printmatrix(std::cout, L3_, "L3_", "--");
+ // L3_.mv(G_tmp,out);
+ // }
+
+ switch (phase)
{
- // printmatrix(std::cout, L2_, "L2_", "--");
- L2_.mv(G_tmp,out);
- }
- else
- {
- // printmatrix(std::cout, L3_, "L3_", "--");
- L3_.mv(G_tmp,out);
+ case 1:
+ {
+ // printmatrix(std::cout, L1_, "L1_", "--");
+ L1_.mv(G_tmp,out);
+ break;
+ }
+ case 2:
+ {
+ // printmatrix(std::cout, L2_, "L2_", "--");
+ L2_.mv(G_tmp,out);
+ break;
+ }
+ case 3:
+ {
+ // printmatrix(std::cout, L3_, "L3_", "--");
+ L3_.mv(G_tmp,out);
+ break;
+ }
+ case 4:
+ {
+ // printmatrix(std::cout, L3_, "L3_", "--");
+ L4_.mv(G_tmp,out);
+ break;
+ }
+ default:
+ DUNE_THROW(Dune::Exception, "Phase number not implemented.");
+ break;
}
// printvector(std::cout, out, "out", "--");
- return VoigtToMatrix(out);
+ return VoigtToMatrix(out); // Muss hier stress zurücktransformieren!
+ // return VoigtToStress(out);
// return image as Matrix again
// return {{out[0], (1.0/sqrt(2.0))*out[5], (1.0/sqrt(2.0))*out[4]}, {(1.0/sqrt(2.0))*out[5], out[1], (1.0/sqrt(2.0))*out[3]}, {(1.0/sqrt(2.0))*out[4], (1.0/sqrt(2.0))*out[3], out[2]}};
// return {{out[0], (1.0/2.0)*out[5], (1.0/2.0)*out[4]}, {(1.0/2.0)*out[5], out[1], (1.0/2.0)*out[3]}, {(1.0/2.0)*out[4], (1.0/2.0)*out[3], out[2]}};
@@ -632,15 +701,32 @@ Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Mo
////////////////////////////////////////////////////////////////////////////////////////
MatrixRT prestrain(const int& phase,const Domain& x) const
{
- if (phase == 1)
- return prestrain1_(x);
- else if (phase == 2)
- return prestrain2_(x);
- else if (phase ==3)
- return prestrain3_(x);
- else
- DUNE_THROW(Dune::Exception, "Phase number not implemented.");
-
+ switch (phase)
+ {
+ case 1:
+ {
+ return prestrain1_(x);
+ break;
+ }
+ case 2:
+ {
+ return prestrain2_(x);
+ break;
+ }
+ case 3:
+ {
+ return prestrain3_(x);
+ break;
+ }
+ case 4:
+ {
+ return prestrain3_(x);
+ break;
+ }
+ default:
+ DUNE_THROW(Dune::Exception, "Phase number not implemented.");
+ break;
+ }
}
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 8ade7d03..e1ba37f5 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -9,7 +9,7 @@
#############################################
### Remove/Comment this when running via Python-Script:
#outputPath=/home/stefan/DUNE/dune-microstructure/outputs
-outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
+outputPath=/home/klaus/Desktop/Dune_release/dune-microstructure/outputs
# --- DEBUG (Output) Option:
@@ -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= 3 3 # computes all levels from first to second entry
#############################################
@@ -35,12 +35,13 @@ numLevels= 2 3 # 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 = "parametrized_laminate"
# --- Choose scale ratio gamma:
@@ -50,7 +51,7 @@ gamma=1.0
#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/materials
-geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/materials
+geometryFunctionPath = /home/klaus/Desktop/Dune_release/dune-microstructure/materials
@@ -69,7 +70,7 @@ geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/mate
#############################################
# Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
#############################################
-Solvertype = 2 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
Solver_verbosity = 0 #(default = 2) degree of information for solver output
diff --git a/materials/parametrized_laminate.py b/materials/parametrized_laminate.py
new file mode 100644
index 00000000..c05161f3
--- /dev/null
+++ b/materials/parametrized_laminate.py
@@ -0,0 +1,68 @@
+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 (abs(x[0]) < (theta/2) and x[2] < 0 ):
+ return 1 #Phase1
+ elif (abs(x[0]) > (theta/2) and x[2] > 0 ):
+ return 2 #Phase2
+ elif (abs(x[0]) < (theta/2) and x[2] > 0 ):
+ return 3 #Phase3
+ else :
+ return 4 #Phase4
+
+
+########### 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] # 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
+######################################################################
+
+# --- Number of material phases
+Phases=4
+
+#--- Define different material phases:
+
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [2.0, 0]
+
+#- PHASE 2
+phase2_type="isotropic"
+materialParameters_phase2 = [1.0, 0]
+
+#- PHASE 3
+phase3_type="isotropic"
+materialParameters_phase3 = [2.0, 0]
+
+#- PHASE 4
+phase4_type="isotropic"
+materialParameters_phase4 = [1.0, 0]
+
+
+#--- define prestrain function for each phase
+# (also works with non-constant values)
+def prestrain_phase1(x):
+ return [[2, 0, 0], [0,2,0], [0,0,2]]
+
+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]]
+
+def prestrain_phase4(x):
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
+
--
GitLab