diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 1fbb2ede0009cb59b1d6f229c4e5626ff6834da8..3e9b4e47b8afa565d5fe714ebab7d096dee9d3ed 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -110,7 +110,7 @@ public:
phiOffset_(basis.size())
{
- assemble();
+ // assemble();
//write VTK call here...
}
@@ -132,7 +132,7 @@ public:
///////////////////////////////
- // getter
+ // Getter
///////////////////////////////
const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}
@@ -1204,7 +1204,8 @@ public:
// --- Write Correctos to VTK:
void writeCorrectorsVTK(const int level)
{
- std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/CellProblem-result";
+ std::string baseName = parameterSet_.get("baseName", "CellProblem-result");
+ std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/" + baseName;
std::cout << "vtkOutputName:" << vtkOutputName << std::endl;
Dune::VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index deae82234e1553f7c3c0646f9e11b854e3b4fe01..067b4876c2c9552335dfd21b8cd5d47210c380bb 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -99,10 +99,9 @@ public:
using MatrixPhase = std::function< MatrixRT(const int&) >;
protected:
-
const GridView& gridView_;
- const Dune::ParameterTree& parameterSet_;
- std::string materialFunctionName_;
+ const Dune::ParameterTree parameterSet_;
+ // std::string materialFunctionName_;
// Elasticity tensors (in Voigt notation)
std::vector<VoigtMatrix<double,3> > L_;
@@ -111,11 +110,6 @@ protected:
Python::Module module_;
- // const FuncScalar mu_;
- // const FuncScalar lambda_;
- // const int phases_;
- // Func2Tensor elasticityTensor_;
- // Func2TensorParam elasticityTensor_;
Func2TensorPhase elasticityTensor_;
@@ -146,103 +140,19 @@ public:
///////////////////////////////
// constructor
///////////////////////////////
+ // prestrainedMaterial(const GridView gridView,
+ // const Dune::ParameterTree& parameterSet)
+ // : gridView_(gridView),
+ // parameterSet_(parameterSet)
+ // {
prestrainedMaterial(const GridView gridView,
- const Dune::ParameterTree& parameterSet)
+ const Dune::ParameterTree& parameterSet,
+ Python::Module pyModule)
: gridView_(gridView),
- parameterSet_(parameterSet)
+ parameterSet_(parameterSet),
+ module_(pyModule)
{
- std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
- std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
-
- // setupMaterial(materialFunctionName_); //old-Version
-
- module_ = Python::import(materialFunctionName_); //TODO use this instead?
-
-
-
- // const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
- // const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
-
-
- //elastic moduli in the order (E1,E2,E3,G12,G23,G31,nu12,nu13,nu23)
-
- //-- Walnut see [Dimwoodie; Timber its nature and behavior p.109]
- // const FieldVector<double,9> walnut_parameters={11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37};
-
- //-- Norway spruce see [Dimwoodie; Timber its nature and behavior p.109]
- // const FieldVector<double,9> Nspruce_parameters={10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31};
-
-
-
-
- // L_[0] = orthotropic(e1,e2,e3,walnut_parameters);
- // L_[1] = orthotropic(e1,e2,e3,Nspruce_parameters);
- // L_[2] = orthotropic(e1,e2,e3,Nspruce_parameters);
-
- setup(materialFunctionName_);
- // setup("material");
-
-
- indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
- localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
-
-
- // L_[0] = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
- // L_[1] = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63, 0.37});
- // L_[2] = transversely_isotropic(e1,e2,e3, {10.7e3,710 ,500,0.51 ,0.31});
-
- // L_[0] = isotropic(mu[0],lambda[0]);
- // L_[1] = isotropic(mu[1],lambda[1]);
- // L_[2] = isotropic(mu[2],lambda[2]);
-
-
-
-
- // L_[0] = {{lambda[0]+2.0*mu[0], lambda[0], lambda[0], 0.0 , 0.0, 0.0},
- // {lambda[0], lambda[0]+2*mu[0], lambda[0], 0.0 ,0.0 ,0.0},
- // {lambda[0], lambda[0], lambda[0]+2.0*mu[0], 0.0, 0.0, 0.0},
- // { 0.0 ,0.0, 0.0, 2.0*mu[0], 0.0, 0.0},
- // { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[0], 0.0},
- // { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[0]}
- // };
-
- // L_[1] = {{lambda[1]+2.0*mu[1], lambda[1], lambda[1], 0.0 , 0.0, 0.0},
- // {lambda[1], lambda[1]+2*mu[1], lambda[1], 0.0 ,0.0 ,0.0},
- // {lambda[1], lambda[1], lambda[1]+2.0*mu[1], 0.0, 0.0, 0.0},
- // { 0.0 ,0.0, 0.0, 2.0*mu[1], 0.0, 0.0},
- // { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[1], 0.0},
- // { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[1]}
- // };
-
- // L_[2] = {{lambda[2]+2.0*mu[2], lambda[2], lambda[2], 0.0 , 0.0, 0.0},
- // {lambda[2], lambda[2]+2.0*mu[2], lambda[2], 0.0 ,0.0 ,0.0},
- // {lambda[2], lambda[2], lambda[2]+2.0*mu[2], 0.0, 0.0, 0.0},
- // { 0.0 ,0.0, 0.0, 2.0*mu[2], 0.0, 0.0},
- // { 0.0 ,0.0, 0.0, 0.0, 2.0*mu[2], 0.0},
- // { 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, L_[0], "L_[0]", "--");
- // L_[0].leftmultiply(D_inv);
- // printmatrix(std::cout, L_[0], "L_[0].leftmultiply(D_inv)", "--");
- // L_[0].rightmultiply(D_);
- // printmatrix(std::cout, L_[0], "L_[0].rightmultiply(D_)", "--");
-
-
- // Python::Module module = Python::import(materialFunctionName_);
- // elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
- // elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
- // auto indicatorFunction = Python::make_function<int>(module.get("indicatorFunction"));
- // indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
- // indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
- // module.get("Phases").toC<int>(Phases_);
- // L_[0] = Python::make_function<MatrixRT>(module.get("L1"));
- // L_[1] = Python::make_function<MatrixRT>(module.get("L2"));
- // L_[2] = Python::make_function<MatrixRT>(module.get("L3"));
- // bool isotropic_ = true; // read from module File
- // auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
+ setup();
}
/**
@@ -250,20 +160,19 @@ public:
* use a constant prestrain.
* - Add switch (constant/non-constant) prestrain?
*
- * @param module python module used
* @param phase material phase
* @return Func2Tensor
*/
-Func2Tensor setupPrestrainPhase(Python::Module module, int phase){
+Func2Tensor setupPrestrainPhase(const int phase){
- auto prestrain = Python::make_function<MatrixRT>(module.get("prestrain_phase" + std::to_string(phase)));
+ auto prestrain = Python::make_function<MatrixRT>(module_.get("prestrain_phase" + std::to_string(phase)));
int axis = 0;
double angle = 0;
try
{
- module.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
- module.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
+ module_.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
+ module_.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
}
catch(Dune::Exception&)
{
@@ -285,17 +194,17 @@ Func2Tensor setupPrestrainPhase(Python::Module module, int phase){
// Dune::FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
-Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
+Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
{
std::string phaseType;
- module.get("phase" + std::to_string(phase) + "_type").toC<std::string>(phaseType);
+ module_.get("phase" + std::to_string(phase) + "_type").toC<std::string>(phaseType);
std::cout << "Setup phase "<< phase << " as " << phaseType << " material." << std::endl;
if(phaseType == "isotropic")
{
Dune::FieldVector<double,2> materialParameters(0);
- module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,2>>(materialParameters);
+ module_.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,2>>(materialParameters);
std::cout << "Lame-Parameters (mu,lambda): (" << materialParameters[0] << "," << materialParameters[1] << ") " << std::endl;
return isotropic(materialParameters[0],materialParameters[1]);
}
@@ -326,8 +235,8 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
// printvector(std::cout, frameVector1, "frameVector1: ", "--");
// printvector(std::cout, frameVector2, "frameVector2: ", "--");
// printvector(std::cout, frameVector3, "frameVector3: ", "--");
- module.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
- module.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
+ module_.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
+ module_.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
}
catch(Dune::Exception& e)
{
@@ -339,21 +248,21 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
if(phaseType == "orthotropic")
{
Dune::FieldVector<double,9> materialParameters(0);
- module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,9>>(materialParameters);
+ module_.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,9>>(materialParameters);
// return orthotropic(frameVector1,frameVector2,frameVector3,materialParameters);
return orthotropic(axis,angle,materialParameters);
}
if(phaseType == "transversely_isotropic")
{
Dune::FieldVector<double,5> materialParameters(0);
- module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,5>>(materialParameters);
+ module_.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,5>>(materialParameters);
// return transversely_isotropic(frameVector1,frameVector2,frameVector3,materialParameters);
return transversely_isotropic(axis,angle,materialParameters);
}
if(phaseType == "general_anisotropic")
{
Dune::FieldMatrix<double,6,6> materialParameters(0); //compliance matric
- module.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldMatrix<double,6,6>>(materialParameters);
+ module_.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldMatrix<double,6,6>>(materialParameters);
// return general_anisotropic(frameVector1,frameVector2,frameVector3,materialParameters);
return general_anisotropic(axis,angle,materialParameters);
}
@@ -367,15 +276,15 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
* @brief Setup method that determines the correct stiffness-matrix (constitutive law)
* for each material phase separately.
*
- * @param name name of the material file being used.
*/
- void setup(const std::string name) // can't use materialFunctionName_ here?
+ void setup()
{
- Python::Module module = Python::import(name);
- indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
+ indicatorFunction_ = Python::make_function<int>(module_.get("indicatorFunction"));
+ indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
+ localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
int phases;
- module.get("Phases").toC<int>(phases);
+ module_.get("Phases").toC<int>(phases);
std::cout << "---- Starting material setup... (Number of Phases: "<< phases << ") " << std::endl;
L_.resize(phases);
@@ -383,13 +292,34 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
for (int i=0; i<phases; i++)
{
- L_[i] = setupPhase(module,i+1);
- prestrain_[i] = setupPrestrainPhase(module, i+1);
+ L_[i] = setupPhase(i+1);
+ prestrain_[i] = setupPrestrainPhase(i+1);
}
std::cout << "Material setup done." << std::endl;
}
+ // void setup(const std::string name) // can't use materialFunctionName_ here?
+ // {
+ // Python::Module module = Python::import(name);
+ // indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
+
+ // int phases;
+ // module.get("Phases").toC<int>(phases);
+ // std::cout << "---- Starting material setup... (Number of Phases: "<< phases << ") " << std::endl;
+
+ // L_.resize(phases);
+ // prestrain_.resize(phases);
+
+ // for (int i=0; i<phases; i++)
+ // {
+ // L_[i] = setupPhase(module,i+1);
+ // prestrain_[i] = setupPrestrainPhase(module, i+1);
+ // }
+
+ // std::cout << "Material setup done." << std::endl;
+ // }
+
////////////////////////////////////////////////////////
// MATERIAL CLASSES DEFINITIONS:
@@ -690,9 +620,10 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
MaterialVtkWriter.addVertexData(
indicatorFunction_P1,
Dune::VTK::FieldInfo("indicatorFunction_P1", Dune::VTK::FieldInfo::Type::scalar, 1));
-
- MaterialVtkWriter.write(outputPath + "/MaterialFunctions-level"+ std::to_string(level) );
- std::cout << "wrote data to file:" + outputPath +"/MaterialFunctions-level" + std::to_string(level) << std::endl;
+
+ std::string baseName = parameterSet_.get("baseName", "CellProblem-result");
+ MaterialVtkWriter.write(outputPath + "/" + baseName + "_MaterialFunction-level"+ std::to_string(level) );
+ std::cout << "wrote data to file:" + outputPath + "/" + baseName + "_MaterialFunction-level"+ std::to_string(level) << std::endl;
return;
};
diff --git a/inputs/material_neukamm.py b/inputs/material_neukamm.py
new file mode 100644
index 0000000000000000000000000000000000000000..cb1724b2a476e7a5c19b957ca4bb828c87e52669
--- /dev/null
+++ b/inputs/material_neukamm.py
@@ -0,0 +1,131 @@
+import math
+
+class ParameterSet(dict):
+ def __init__(self, *args, **kwargs):
+ super(ParameterSet, self).__init__(*args, **kwargs)
+ self.__dict__ = self
+
+parameterSet = ParameterSet()
+
+#############################################
+# Paths
+#############################################
+parameterSet.outputPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/outputs'
+parameterSet.baseName= 'material_neukamm' #(needed for Output-Filename)
+
+#############################################
+# Material Setup
+#############################################
+
+# --- Choose scale ratio gamma:
+parameterSet.gamma = 1.0
+
+# --- Number of material phases
+Phases = 3
+
+#--- Indicator function for material phases
+# x[0] : y1-component
+# x[1] : y2-component
+# x[2] : x3-component
+#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
+
+########### 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
+######################################################################
+
+#--- 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]
+
+#--- 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]]
+
+
+
+#############################################
+# Grid parameters
+#############################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+parameterSet.numLevels= '3 3' # computes all levels from first to second entry
+
+#############################################
+# Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1 #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
+#parameterSet.arbitraryLocalIndex = 7 #(default = 0)
+#parameterSet.arbitraryElementNumber = 3 #(default = 0)
+
+#############################################
+# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0 #(default = 2) degree of information for solver output
+
+
+#############################################
+# Write/Output options #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1 # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1 # VTK norm of B (currently not implemented)
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0 #(default=false)
+
+# --- Write Correctos to VTK-File:
+parameterSet.write_VTK = 1
+
+# --- (Optional output) L2Error, integral mean:
+#parameterSet.write_L2Error = 1
+#parameterSet.write_IntegralMean = 1
+
+# --- check orthogonality (75) from paper:
+parameterSet.write_checkOrthogonality = 1
+
+# --- Write corrector-coefficients to log-File:
+#parameterSet.write_corrector_phi1 = 1
+#parameterSet.write_corrector_phi2 = 1
+#parameterSet.write_corrector_phi3 = 1
+
+# --- Print Condition number of matrix (can be expensive):
+#parameterSet.print_conditionNumber= 1 #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+#parameterSet.write_toMATLAB = 0 # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/inputs/material_orthotropic.py b/inputs/material_orthotropic.py
new file mode 100644
index 0000000000000000000000000000000000000000..4c140c6e26fdaeb58611ef47ff26a47c7a115443
--- /dev/null
+++ b/inputs/material_orthotropic.py
@@ -0,0 +1,130 @@
+import math
+import numpy as np
+
+class ParameterSet(dict):
+ def __init__(self, *args, **kwargs):
+ super(ParameterSet, self).__init__(*args, **kwargs)
+ self.__dict__ = self
+
+parameterSet = ParameterSet()
+
+#############################################
+# Paths
+#############################################
+parameterSet.outputPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/outputs'
+parameterSet.baseName= 'material_orthotropic' #(needed for Output-Filename)
+
+#############################################
+# Material Setup
+#############################################
+
+# --- Choose scale ratio gamma:
+parameterSet.gamma = 1.0
+
+# --- Number of material phases
+Phases = 3
+
+#--- Indicator function for material phases
+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] # 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
+######################################################################
+
+#--- 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]
+materialParameters_phase2 = [11.2e3,630,1190,700,230,960,0.63 ,0.49,0.37]
+
+phase2_axis = 2
+phase2_angle = np.pi/2.0
+
+#- PHASE 3
+phase3_type="isotropic"
+materialParameters_phase3 = [60, 25]
+
+
+#--- 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]]
+
+
+
+#############################################
+# Grid parameters
+#############################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+parameterSet.numLevels= '3 3' # computes all levels from first to second entry
+
+#############################################
+# Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1 #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
+#parameterSet.arbitraryLocalIndex = 7 #(default = 0)
+#parameterSet.arbitraryElementNumber = 3 #(default = 0)
+
+#############################################
+# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0 #(default = 2) degree of information for solver output
+
+
+#############################################
+# Write/Output options #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1 # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1 # VTK norm of B (currently not implemented)
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0 #(default=false)
+
+# --- Write Correctos to VTK-File:
+parameterSet.write_VTK = 1
+
+# --- (Optional output) L2Error, integral mean:
+#parameterSet.write_L2Error = 1
+#parameterSet.write_IntegralMean = 1
+
+# --- check orthogonality (75) from paper:
+parameterSet.write_checkOrthogonality = 1
+
+# --- Write corrector-coefficients to log-File:
+#parameterSet.write_corrector_phi1 = 1
+#parameterSet.write_corrector_phi2 = 1
+#parameterSet.write_corrector_phi3 = 1
+
+# --- Print Condition number of matrix (can be expensive):
+#parameterSet.print_conditionNumber= 1 #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+#parameterSet.write_toMATLAB = 0 # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/inputs/parametrized_laminate.py b/inputs/parametrized_laminate.py
new file mode 100644
index 0000000000000000000000000000000000000000..0ca84a2a83a689ddbe6f3edb97e0a1f02f8843eb
--- /dev/null
+++ b/inputs/parametrized_laminate.py
@@ -0,0 +1,133 @@
+import math
+
+class ParameterSet(dict):
+ def __init__(self, *args, **kwargs):
+ super(ParameterSet, self).__init__(*args, **kwargs)
+ self.__dict__ = self
+
+parameterSet = ParameterSet()
+
+#############################################
+# Paths
+#############################################
+parameterSet.outputPath = '/home/klaus/Desktop/Dune_release/dune-microstructure/outputs'
+parameterSet.baseName= 'parametrized_laminate' #(needed for Output-Filename)
+
+#############################################
+# Material Setup
+#############################################
+
+# --- Choose scale ratio gamma:
+parameterSet.gamma = 1.0
+
+# --- Number of material phases
+Phases = 4
+
+#--- Indicator function for material phases
+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
+######################################################################
+
+#--- 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]]
+
+
+
+#############################################
+# Grid parameters
+#############################################
+## numLevels : Number of Levels on which solution is computed. starting with a 2x2x2 cube mesh.
+## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
+#----------------------------------------------------
+parameterSet.numLevels= '3 3' # computes all levels from first to second entry
+
+#############################################
+# Assembly options
+#############################################
+parameterSet.set_IntegralZero = 1 #(default = false)
+parameterSet.set_oneBasisFunction_Zero = 1 #(default = false)
+#parameterSet.arbitraryLocalIndex = 7 #(default = 0)
+#parameterSet.arbitraryElementNumber = 3 #(default = 0)
+
+#############################################
+# Solver Options, Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
+#############################################
+parameterSet.Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+parameterSet.Solver_verbosity = 0 #(default = 2) degree of information for solver output
+
+
+#############################################
+# Write/Output options #(default=false)
+#############################################
+# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
+parameterSet.write_materialFunctions = 1 # VTK indicator function for material/prestrain definition
+#parameterSet.write_prestrainFunctions = 1 # VTK norm of B (currently not implemented)
+
+# --- (Additional debug output)
+parameterSet.print_debug = 0 #(default=false)
+
+# --- Write Correctos to VTK-File:
+parameterSet.write_VTK = 1
+
+# --- (Optional output) L2Error, integral mean:
+#parameterSet.write_L2Error = 1
+#parameterSet.write_IntegralMean = 1
+
+# --- check orthogonality (75) from paper:
+parameterSet.write_checkOrthogonality = 1
+
+# --- Write corrector-coefficients to log-File:
+#parameterSet.write_corrector_phi1 = 1
+#parameterSet.write_corrector_phi2 = 1
+#parameterSet.write_corrector_phi3 = 1
+
+# --- Print Condition number of matrix (can be expensive):
+#parameterSet.print_conditionNumber= 1 #(default=false)
+
+# --- write effective quantities to Matlab-folder for symbolic minimization:
+#parameterSet.write_toMATLAB = 0 # writes effective quantities to .txt-files QMatrix.txt and BMatrix.txt
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index ea5d8668411dac41e4cf8894f95c64f893f80a18..0877756bbff3dcd00fe9648424c3c02c67bfef88 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -3,6 +3,8 @@
#include <vector>
#include <fstream>
+#include <fenv.h>
+
#include <iostream>
#include <dune/common/indices.hh>
#include <dune/common/bitsetvector.hh>
@@ -116,51 +118,69 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char *argv[])
{
- MPIHelper::instance(argc, argv);
-
- Dune::Timer globalTimer;
+ feenableexcept(FE_INVALID);
+ MPIHelper::instance(argc, argv);
+
+ Dune::Timer globalTimer;
- ParameterTree parameterSet;
- if (argc < 2)
- ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
- else
- {
- ParameterTreeParser::readINITree(argv[1], parameterSet);
+
+ if (argc < 3)
+ DUNE_THROW(Exception, "Usage: ./Cell-Problem <python path> <python module without extension>");
+
+ // Start Python interpreter
+ Python::start();
+ auto pyMain = Python::main();
+ pyMain.runStream()
+ << std::endl << "import math"
+ << std::endl << "import sys"
+ << std::endl << "sys.path.append('" << argv[1] << "')" << std::endl;
+ auto pyModule = pyMain.import(argv[2]);
+
+ ParameterTree parameterSet;
+ pyModule.get("parameterSet").toC(parameterSet);
+ // read possible further parameters from the command line
ParameterTreeParser::readOptions(argc, argv, parameterSet);
- }
+ // Print all parameters, to make them appear in the log file
+ std::cout << "Executable: harmonicmaps, with parameters:" << std::endl;
+ parameterSet.report();
+
+// ParameterTree parameterSet;
+// if (argc < 2)
+// ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
+// else
+// {
+// ParameterTreeParser::readINITree(argv[1], parameterSet);
+// ParameterTreeParser::readOptions(argc, argv, parameterSet);
+// }
//--- Output setter
- std::string outputPath = parameterSet.get("outputPath", "../../outputs");
-
+ std::string baseName = parameterSet.get("baseName", "CellProblem-result");
+ std::string outputPath = parameterSet.get("outputPath", "../../outputs") ;
//--- setup Log-File
std::fstream log;
- log.open(outputPath + "/output.txt" ,std::ios::out);
+ log.open(outputPath + "/" + baseName + "_log.txt" ,std::ios::out);
- std::cout << "outputPath:" << outputPath << std::endl;
// parameterSet.report(log); // short Alternativ
//--- Get Path for Material/Geometry functions
- std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath");
- //--- Start Python interpreter
- Python::start();
- Python::Reference main = Python::import("__main__");
- Python::run("import math");
- Python::runStream()
- << std::endl << "import sys"
- << std::endl << "sys.path.append('" << geometryFunctionPath << "')"
- << std::endl;
+ // std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath");
+ // //--- Start Python interpreter
+ // Python::start();
+ // Python::Reference main = Python::import("__main__");
+ // Python::run("import math");
+ // Python::runStream()
+ // << std::endl << "import sys"
+ // << std::endl << "sys.path.append('" << geometryFunctionPath << "')"
+ // << std::endl;
constexpr int dim = 3;
- constexpr int dimWorld = 3;
// Debug/Print Options
bool print_debug = parameterSet.get<bool>("print_debug", false);
-
// VTK-write options
- bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
- bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false);
+// bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false); //Not implemented yet.
///////////////////////////////////
// Generate the grid
@@ -169,16 +189,17 @@ int main(int argc, char *argv[])
FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
- std::array<int,2> numLevels = parameterSet.get<std::array<int,2>>("numLevels", {1,3});
+ std::array<int,2> numLevels = parameterSet.get<std::array<int,2>>("numLevels", {3,3});
int levelCounter = 0;
///////////////////////////////////
- // Create Date Storage
+ // Create Data Storage
///////////////////////////////////
//--- Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
- std::vector<std::variant<std::string, size_t , double>> Storage_Error;
- //--- Storage:: | level | q1 | q2 | q3 | q12 | q23 | b1 | b2 | b3 |
+// std::vector<std::variant<std::string, size_t , double>> Storage_Error;
+
+ //--- Storage:: | level | q1 | q2 | q3 | q12 | q13 | q23 | b1 | b2 | b3 |
std::vector<std::variant<std::string, size_t , double>> Storage_Quantities;
@@ -189,7 +210,7 @@ int main(int argc, char *argv[])
std::cout << "GridLevel: " << level << std::endl;
std::cout << " ----------------------------------" << std::endl;
- Storage_Error.push_back(level);
+ // Storage_Error.push_back(level);
Storage_Quantities.push_back(level);
std::array<int, dim> nElements = {(int)std::pow(2,level) ,(int)std::pow(2,level) ,(int)std::pow(2,level)};
std::cout << "Number of Grid-Elements in each direction: " << nElements << std::endl;
@@ -206,7 +227,7 @@ int main(int argc, char *argv[])
using namespace Functions::BasisFactory;
Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
- //--- get PeriodicIndices for periodicBasis (Don't do the following in real life: It has quadratic run-time in the number of vertices.)
+ //--- Get PeriodicIndices for periodicBasis (Don't do the following in real life: It has quadratic run-time in the number of vertices.)
for (const auto& v1 : vertices(gridView_CE))
for (const auto& v2 : vertices(gridView_CE))
if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
@@ -217,7 +238,7 @@ int main(int argc, char *argv[])
//--- setup first order periodic Lagrange-Basis
auto Basis_CE = makeBasis(
gridView_CE,
- power<dim>( // eig dimworld?!?!
+ power<dim>(
Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
flatLexicographic()
//blockedInterleaved() // Not Implemented
@@ -229,7 +250,7 @@ int main(int argc, char *argv[])
///////////////////////////////////
// Create prestrained material object
///////////////////////////////////
- auto material_ = prestrainedMaterial(gridView_CE,parameterSet);
+ auto material = prestrainedMaterial(gridView_CE,parameterSet,pyModule);
// --- get scale ratio
double gamma = parameterSet.get<double>("gamma",1.0);
@@ -237,9 +258,8 @@ int main(int argc, char *argv[])
//------------------------------------------------------------------------------------------------
//--- Compute Correctors
- // auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
- // auto correctorComputer = CorrectorComputer(Basis_CE, material_, muTerm, lambdaTerm, gamma, log, parameterSet);
- auto correctorComputer = CorrectorComputer(Basis_CE, material_, gamma, log, parameterSet);
+ auto correctorComputer = CorrectorComputer(Basis_CE, material, gamma, log, parameterSet);
+ correctorComputer.assemble();
correctorComputer.solve();
//--- Check Correctors (options):
@@ -255,13 +275,13 @@ int main(int argc, char *argv[])
correctorComputer.checkSymmetry();
//--- Compute effective quantities
- auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material_);
+ auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material);
effectiveQuantitiesComputer.computeEffectiveQuantities();
//--- write material indicator function to VTK
- if (write_materialFunctions)
+ if (parameterSet.get<bool>("write_materialFunctions", false))
{
- material_.writeVTKMaterialFunctions(nElements,level);
+ material.writeVTKMaterialFunctions(nElements,level);
}
//--- TEST:: Compute Qeff without using the orthogonality (75)...
@@ -291,6 +311,7 @@ int main(int argc, char *argv[])
Storage_Quantities.push_back(Qeff[1][1] );
Storage_Quantities.push_back(Qeff[2][2] );
Storage_Quantities.push_back(Qeff[0][1] );
+ Storage_Quantities.push_back(Qeff[0][2] );
Storage_Quantities.push_back(Qeff[1][2] );
Storage_Quantities.push_back(Beff[0]);
Storage_Quantities.push_back(Beff[1]);
@@ -301,6 +322,7 @@ int main(int argc, char *argv[])
log << "q2=" << Qeff[1][1] << std::endl;
log << "q3=" << Qeff[2][2] << std::endl;
log << "q12=" << Qeff[0][1] << std::endl;
+ log << "q13=" << Qeff[0][2] << std::endl;
log << "q23=" << Qeff[1][2] << std::endl;
log << std::fixed << std::setprecision(6) << "q_onetwo=" << Qeff[0][1] << std::endl;
log << "b1=" << Beff[0] << std::endl;
@@ -321,22 +343,24 @@ int main(int argc, char *argv[])
<< center("q2",tableWidth) << " | "
<< center("q3",tableWidth) << " | "
<< center("q12",tableWidth) << " | "
+ << center("q13",tableWidth) << " | "
<< center("q23",tableWidth) << " | "
<< center("b1",tableWidth) << " | "
<< center("b2",tableWidth) << " | "
<< center("b3",tableWidth) << " | " << "\n";
- std::cout << std::string(tableWidth*9 + 3*9, '-') << "\n";
- log << std::string(tableWidth*9 + 3*9, '-') << "\n";
+ std::cout << std::string(tableWidth*10 + 3*10, '-') << "\n";
+ log << std::string(tableWidth*10 + 3*10, '-') << "\n";
log << center("Levels ",tableWidth) << " | "
<< center("q1",tableWidth) << " | "
<< center("q2",tableWidth) << " | "
<< center("q3",tableWidth) << " | "
<< center("q12",tableWidth) << " | "
+ << center("q13",tableWidth) << " | "
<< center("q23",tableWidth) << " | "
<< center("b1",tableWidth) << " | "
<< center("b2",tableWidth) << " | "
<< center("b3",tableWidth) << " | " << "\n";
- log << std::string(tableWidth*9 + 3*9, '-') << "\n";
+ log << std::string(tableWidth*10 + 3*10, '-') << "\n";
int StorageCount2 = 0;
for(auto& v: Storage_Quantities)
@@ -344,16 +368,16 @@ int main(int argc, char *argv[])
std::visit([tableWidth](auto&& arg){std::cout << center(prd(arg,5,1),tableWidth) << " | ";}, v);
std::visit([tableWidth, &log](auto&& arg){log << center(prd(arg,5,1),tableWidth) << " & ";}, v);
StorageCount2++;
- if(StorageCount2 % 9 == 0 )
+ if(StorageCount2 % 10 == 0 )
{
std::cout << std::endl;
log << std::endl;
}
}
- std::cout << std::string(tableWidth*9 + 3*9, '-') << "\n";
- log << std::string(tableWidth*9 + 3*9, '-') << "\n";
+ std::cout << std::string(tableWidth*10 + 3*10, '-') << "\n";
+ log << std::string(tableWidth*10 + 3*10, '-') << "\n";
- log.close(); //close log-file
+ log.close();
std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
}