From 946583ec766b5cd3b831925f1dffb1831c80233c Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Wed, 31 Aug 2022 19:49:44 +0200
Subject: [PATCH] adjust rest of code
---
dune/microstructure/CorrectorComputer.hh | 23 +-
.../EffectiveQuantitiesComputer.hh | 47 ++--
dune/microstructure/prestrainedMaterial.hh | 249 +++++++++---------
geometries/material.py | 61 +++--
inputs/cellsolver.parset | 7 +-
5 files changed, 198 insertions(+), 189 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index cad9d1ae..93c6973d 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -255,7 +255,7 @@ public:
elementMatrix = 0;
- auto elasticityTensor = material_.getElasticityTensor();
+ // auto elasticityTensor = material_.getElasticityTensor();
// auto indicatorFunction = material_.getIndicatorFunction();
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
localIndicatorFunction.bind(element);
@@ -375,11 +375,11 @@ public:
// std::cout << "--------------- \n";
// printmatrix(std::cout, defGradientU, "defGradientU", "--");
- // printmatrix(std::cout, material_.applyL(defGradientU,localIndicatorFunction(quadPos)), "material_.applyL(defGradientU,localIndicatorFunction(quadPos))", "--");
+ // printmatrix(std::cout, material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos))", "--");
// printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
- double energyDensity= scalarProduct(material_.applyL(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
+ double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
@@ -420,7 +420,7 @@ public:
- double energyDensityGphi = scalarProduct(material_.applyL(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
+ double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) <<std::endl;
@@ -469,7 +469,7 @@ public:
- double energyDensityGG = scalarProduct(material_.applyL(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
+ double energyDensityGG = scalarProduct(material_.applyElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
@@ -522,7 +522,7 @@ public:
// using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
- auto elasticityTensor = material_.getElasticityTensor();
+ // auto elasticityTensor = material_.getElasticityTensor();
// auto indicatorFunction = material_.getIndicatorFunction();
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
localIndicatorFunction.bind(element);
@@ -607,7 +607,7 @@ public:
- double energyDensity= scalarProduct(material_.applyL((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
+ double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
@@ -633,7 +633,7 @@ public:
// double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
- double energyDensityfG= scalarProduct(material_.applyL((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
+ double energyDensityfG= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
@@ -1082,8 +1082,8 @@ public:
// std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
log_ << "Solver-type used: " <<" UMFPACK-Solver" << std::endl;
}
- std::cout << "Finished solving Corrector problems!" << std::endl;
- std::cout << "Time for solving:" << SolverTimer.elapsed() << std::endl;
+ std::cout << "--- Corrector computation finished ---" << std::endl;
+ std::cout << "Time required for solving:" << SolverTimer.elapsed() << std::endl;
////////////////////////////////////////////////////////////////////////////////////
// Extract phi_alpha & M_alpha coefficients
@@ -1409,7 +1409,8 @@ public:
auto tmp = G + *mContainer[a] + strain1;
- double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
+ double energyDensity= scalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
+ // double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index d704637f..0d56aade 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -39,15 +39,15 @@ public:
protected:
CorrectorComputer<Basis,Material>& correctorComputer_;
// Func2Tensor prestrain_;
- MatrixPhase prestrain_;
+ // MatrixPhase prestrain_;
const Material& material_;
// Func2int indicatorFunction_;
public:
- MatrixRT Qeff_; //effective moduli
+ MatrixRT Qeff_; // Effective moduli
VectorRT Bhat_;
- VectorRT Beff_;
+ VectorRT Beff_; // Effective prestrain
///////////////////////////////
@@ -55,11 +55,13 @@ public:
///////////////////////////////
EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
const Material& material)
- : correctorComputer_(correctorComputer),
- material_(material)
+ : correctorComputer_(correctorComputer),
+ material_(material)
{
- prestrain_ = material_.getPrestrainFunction();
+ // prestrain_ = material_.getPrestrainFunction();
+ std::cout << "starting effective quantities computation..." << std::endl;
+
}
@@ -129,22 +131,22 @@ public:
// using GridView = typename Basis::GridView;
- for (const auto& e : elements(basis.gridView()))
+ for (const auto& element : elements(basis.gridView()))
{
- localView.bind(e);
- matrixFieldG1.bind(e);
- matrixFieldG2.bind(e);
- localfun_a.bind(e);
+ localView.bind(element);
+ matrixFieldG1.bind(element);
+ matrixFieldG2.bind(element);
+ localfun_a.bind(element);
// DerPhi2.bind(e);
// mu.bind(e);
// lambda.bind(e);
// prestrainFunctional.bind(e);
- localIndicatorFunction.bind(e);
+ localIndicatorFunction.bind(element);
double elementEnergy = 0.0;
double elementPrestrain = 0.0;
- auto geometry = e.geometry();
+ auto geometry = element.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
@@ -153,7 +155,7 @@ public:
// int orderQR = 1;
// int orderQR = 2;
// int orderQR = 3;
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quad)
{
@@ -171,13 +173,24 @@ public:
auto X1 = G1 + Chi1;
// auto X2 = G2 + Chi2;
- double energyDensity= scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
+ double energyDensity= scalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
+ // double energyDensity= scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
if (b==0)
{
- auto prestrainPhaseValue = prestrain_(localIndicatorFunction(quadPos));
- elementPrestrain += scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue)) * quadPoint.weight() * integrationElement;
+ // auto prestrainPhaseValue_old = prestrain_(localIndicatorFunction(quadPos));
+ auto prestrainPhaseValue = material_.prestrain(localIndicatorFunction(quadPos),element.geometry().global(quadPos));
+ // if (localIndicatorFunction(quadPos) != 3)
+ // {
+ // std::cout << "element.geometry().global(quadPos):"<< element.geometry().global(quadPos) << std::endl;
+ // printmatrix(std::cout, prestrainPhaseValue_old, "prestrainPhaseValue_old", "--");
+ // printmatrix(std::cout, prestrainPhaseValue, "prestrainPhaseValue", "--");
+ // }
+ auto value = scalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue));
+
+ elementPrestrain += value * quadPoint.weight() * integrationElement;
+ // elementPrestrain += scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue)) * quadPoint.weight() * integrationElement;
// elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos)) * quadPoint.weight() * integrationElement;
}
}
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index e158dac7..3b7c1a57 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -113,6 +113,14 @@ protected:
FieldMatrix<double,6,6> L1_;
FieldMatrix<double,6,6> L2_;
FieldMatrix<double,6,6> L3_;
+ FieldMatrix<double,6,6> L4_; //not used yet
+
+ Func2Tensor prestrain1_;
+ Func2Tensor prestrain2_;
+ Func2Tensor prestrain3_;
+ Func2Tensor prestrain4_; //not used yet
+
+ Python::Module module_;
// const FuncScalar mu_;
// const FuncScalar lambda_;
@@ -129,20 +137,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)}
+ // };
@@ -158,22 +166,23 @@ public:
std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
- setupMaterial(materialFunctionName_);
+ // setupMaterial(materialFunctionName_); //old-Version
- indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
- localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
+ 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};
+
+
+ // 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};
+ // 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};
+ // const FieldVector<double,9> Nspruce_parameters={10.7e3,430,710,620,23,500, 0.51 ,0.38,0.31};
//frame vectors
VectorRT e1 = {1.0,0.0,0.0};
@@ -186,9 +195,12 @@ public:
// L2_ = orthotropic(e1,e2,e3,Nspruce_parameters);
// L3_ = orthotropic(e1,e2,e3,Nspruce_parameters);
+ setup(materialFunctionName_);
+ // setup("material");
- setup("material");
+ indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
+ localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
// L1_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
@@ -251,10 +263,6 @@ public:
- ///////////////////////////////////////////////////////////////////////
-
-// void elasticitySetter //set L1_...
-
FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module module, int phase)
{
@@ -325,94 +333,93 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
}
else
DUNE_THROW(Exception, " Unknown material class for phase ");
-
}
+ //---function that determines elasticity Tensor
+ // void setup(const std::string name, const int Phases) // cant use materialFunctionName_ here!?
+ void setup(const std::string name) // cant use materialFunctionName_ here!?
+ {
- //---function that determines elasticity Tensor
- // void setup(const std::string name, const int Phases) // cant use materialFunctionName_ here!?
- void setup(const std::string name) // cant use materialFunctionName_ here!?
- {
-
- Python::Module module = Python::import(name);
- indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
+ Python::Module module = Python::import(name);
+ indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
- std::cout << "---- setup material... " << std::endl;
-
- int Phases;
- module.get("Phases").toC<int>(Phases);
- std::cout << "Number of Phases: " << Phases << std::endl;
-
-
- switch (Phases)
- {
- case 1:
- {
- std::string phase1_type;
- module.get("phase1_type").toC<std::string>(phase1_type);
- L1_ = setupPhase(phase1_type, module,1);
- break;
- }
- case 2:
- {
- 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::cout << "---- setup material... " << std::endl;
+
+ int Phases;
+ module.get("Phases").toC<int>(Phases);
+ std::cout << "Number of Phases: " << Phases << std::endl;
- L1_ = setupPhase(phase1_type, module,1);
- L2_ = setupPhase(phase2_type, module,2);
- break;
- }
- case 3:
- {
- 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 phase1_type = parameterSet_.get<std::string>("phase1_type", "isotropic");
- // std::string phase2_type = parameterSet_.get<std::string>("phase2_type", "isotropic");
- // std::string phase3_type = parameterSet_.get<std::string>("phase3_type", "isotropic");
-
- L1_ = setupPhase(phase1_type, module,1);
- L2_ = setupPhase(phase2_type, module,2);
- L3_ = setupPhase(phase3_type, module,3);
-
- // 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;
- }
- default:
- break;
- }
- std::cout << "Material setup done." << std::endl;
+ switch (Phases)
+ {
+ case 1:
+ {
+ std::string phase1_type;
+ module.get("phase1_type").toC<std::string>(phase1_type);
+ L1_ = setupPhase(phase1_type, module,1);
+ prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
+ break;
}
+ case 2:
+ {
+ 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);
- ///////////////////////////////////////////////////////////////////////
-
+ L1_ = setupPhase(phase1_type, module,1);
+ L2_ = setupPhase(phase2_type, module,2);
+ prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
+ prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
+ break;
+ }
+ case 3:
+ {
+ 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);
+
+ L1_ = setupPhase(phase1_type, module,1);
+ L2_ = setupPhase(phase2_type, module,2);
+ L3_ = setupPhase(phase3_type, module,3);
+
+ 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"));
+
+ // 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;
+ }
+ default:
+ break;
+ }
+ std::cout << "Material setup done." << std::endl;
+ }
+ ///////////////////////////////////////////////////////////////////////
//---function that determines elasticity Tensor
@@ -586,19 +593,14 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
- MatrixRT applyL(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]};
FieldVector<double,6> G_tmp = MatrixToVoigt(G);
FieldVector<double,6> out(0);
-
// printvector(std::cout, G_tmp, "G_tmp", "--");
- //
- // printmatrix(std::cout, L2_, "L2_", "--");
- // printmatrix(std::cout, L3_, "L3_", "--");
-
if (phase == 1)
{
// printmatrix(std::cout, L1_, "L1_", "--");
@@ -624,28 +626,17 @@ FieldMatrix<double,6,6> setupPhase(const std::string phaseType, Python::Module m
// return {{out[0], out[5], out[4]}, {out[5], out[1], out[3]}, {out[4], out[3], out[2]}};
}
- // MatrixRT prestrain(const int& phase,const Domain& x) const
- // {
-
-
- // if (phase == 1)
- // {
- // // printmatrix(std::cout, L1_, "L1_", "--");
- // prestrain1_.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);
- // }
-
- // }
+ ////////////////////////////////////////////////////////////////////////////////////////
+ MatrixRT prestrain(const int& phase,const Domain& x) const
+ {
+ if (phase == 1)
+ return prestrain1_(x);
+ else if (phase == 2)
+ return prestrain2_(x);
+ else
+ return prestrain3_(x);
+ }
diff --git a/geometries/material.py b/geometries/material.py
index d1cc5ade..b2756eb6 100644
--- a/geometries/material.py
+++ b/geometries/material.py
@@ -3,11 +3,11 @@ from python_matrix_operations import *
import ctypes
import os
import sys
+#---------------------------------------------------------------
-
-
+#--- define indicator function
def indicatorFunction(x):
theta=0.25
factor=1
@@ -22,48 +22,50 @@ def indicatorFunction(x):
########### Options for material phases: #################################
# 1. "isotropic" 2. "orthotropic" 3. "transversely_isotropic" 4. "general_anisotropic"
#########################################################################
-## Notation (materialParameters_phase):
-# 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
+## 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:
-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] # Nspruce_parameters (values for compliance matrix)
+#--- Define different material phases:
-phase2_type="transversely_isotropic"
-materialParameters_phase2 = [11.2e3,1190,960,0.63 ,0.37]
+#- PHASE 1
+phase1_type="isotropic"
+materialParameters_phase1 = [80, 80]
-# phase1_type="isotropic"
-# materialParameters_phase1 = [80, 80]
-#
-# phase2_type="isotropic"
-# materialParameters_phase2 = [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)
-# phase3_type="isotropic"
-# materialParameters_phase3 = [60, 25]
+#- PHASE 2
+# phase2_type="transversely_isotropic"
+# materialParameters_phase2 = [11.2e3,1190,960,0.63 ,0.37]
+phase2_type="isotropic"
+materialParameters_phase2 = [80, 80]
-# 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]])
+#- 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
def prestrain_phase1(x):
return [[1, 0, 0], [0,1,0], [0,0,1]]
@@ -72,3 +74,4 @@ def prestrain_phase2(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/inputs/cellsolver.parset b/inputs/cellsolver.parset
index c88beabb..9ce849f0 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -13,7 +13,7 @@ outputPath=/home/klaus/Desktop/Dune-Testing/dune-microstructure/outputs
# --- DEBUG (Output) Option:
-print_debug = true #(default=false)
+#print_debug = true #(default=false)
@@ -36,7 +36,8 @@ numLevels= 2 3 # computes all levels from first to second entry
# --- Choose material definition:
#materialFunction = "two_phase_material_1"
#materialFunction = "two_phase_material_2"
-materialFunction = "material_neukamm"
+#materialFunction = "material_neukamm"
+materialFunction = "material"
#materialFunction = "material_2"
#materialFunction = "homogeneous"
@@ -66,7 +67,7 @@ geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geom
#############################################
# Solver Type: #1: CG - SOLVER , #2: GMRES - SOLVER, #3: QR - SOLVER (default), #4: UMFPACK - SOLVER
#############################################
-Solvertype = 3 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
+Solvertype = 2 # recommended to use iterative solver (e.g GMRES) for finer grid-levels
Solver_verbosity = 0 #(default = 2) degree of information for solver output
--
GitLab