diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index f4b52c9c0c251b1569db902b9c4b7f8539d06e07..180552f301d6ec515bdd517ee95a24f3baf0c589 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -38,6 +38,7 @@ public:
using FuncScalar = std::function< ScalarRT(const Domain&) >;
using FuncVector = std::function< VectorRT(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ using Func2int = std::function< int(const Domain&) >;
using VectorCT = BlockVector<FieldVector<double,1> >;
using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
@@ -183,6 +184,8 @@ public:
shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
+
+
// Get the occupation pattern of the stiffness matrix
void getOccupationPattern(MatrixIndexSet& nb)
{
@@ -266,10 +269,13 @@ public:
auto elasticityTensor = material_.getElasticityTensor();
- auto indicatorFunction = material_.getIndicatorFunction();
+ // auto indicatorFunction = material_.getIndicatorFunction();
+ // auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ // localIndicatorFunction.bind(element);
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF);
// indicatorFunction.bind(element);
+ // Func2int indicatorFunction = *material_.getIndicatorFunction();
// LocalBasis-Offset
@@ -350,9 +356,34 @@ public:
// double energyDensity= scalarProduct(etmp,defGradientV);
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),defGradientV);
- double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
+ // auto tmp_1 = scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
+ // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+ // if (abs(tmp_1-tmp_2)>1e-2)
+ // {
+ // std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
+ // std::cout << "scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))<< std::endl;
+ // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)<< std::endl;
+ // }
+
+
+ // Python::Module module = Python::import("material_neukamm");
+ // auto PythonindicatorFunction = Python::make_function<int>(module.get("f"));
+ // if (PythonindicatorFunction(element.geometry().global(quadPos)) != material_.applyIndicatorFunction(element.geometry().global(quadPos)))
+ // {
+ // std::cout <<" element.geometry().global(quadPos): " << element.geometry().global(quadPos) << std::endl;
+ // std::cout << "PythonindicatorFunction(element.geometry().global(quadPos)): " << PythonindicatorFunction(element.geometry().global(quadPos)) << std::endl;
+ // // std::cout << "mu(quadPos):" << mu(quadPos) << std::endl;
+ // std::cout << "indicatorFunction(element.geometry().global(quadPos)): " << indicatorFunction(element.geometry().global(quadPos)) << std::endl;
+ // std::cout << "material_.applyIndicatorFunction(element.geometry().global(quadPos)): " << material_.applyIndicatorFunction(element.geometry().global(quadPos)) << std::endl;
+ // std::cout << "indicatorFunction_material_1: " << indicatorFunction_material_1(element.geometry().global(quadPos)) << std::endl;
+ // }
+ double energyDensity= scalarProduct(material_.ElasticityTensorGlobal(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
+ // double energyDensity= scalarProduct(material_.ElasticityTensorPhase(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));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
+ // double energyDensity= scalarProduct(material_.localElasticityTensor(defGradientU,quadPos,element),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal(defGradientU,quadPos),defGradientV);
@@ -369,12 +400,29 @@ public:
// "m*phi" & "phi*m" - part
for( size_t m=0; m<3; m++)
{
-
+
+
+ // auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) ;
+ // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+ // if (abs(tmp_1-tmp_2)>1e-2)
+ // {
+ // std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
+ // 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;
+ // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
+ // }
+
+
+ double energyDensityGphi = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
+ // double energyDensityGphi = scalarProduct(material_.ElasticityTensorPhase(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;
+ // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
- double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
+ // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(defGradientV));
+ // double energyDensityGphi= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),defGradientV);
// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
@@ -399,12 +447,25 @@ public:
// std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
// std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
- double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
+ // auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
+ // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
+ // if (abs(tmp_1-tmp_2)>1e-2)
+ // {
+ // std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
+ // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])):" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])) << std::endl;
+ // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])<< std::endl;
+
+ // }
+ double energyDensityGG = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
+ // double energyDensityGG = scalarProduct(material_.ElasticityTensorPhase(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]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(basisContainer[n]));
+ // double energyDensityGG= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),basisContainer[n]);
// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
@@ -444,11 +505,15 @@ public:
// using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
auto elasticityTensor = material_.getElasticityTensor();
- auto indicatorFunction = material_.getIndicatorFunction();
+ // auto indicatorFunction = material_.getIndicatorFunction();
+ // auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ // localIndicatorFunction.bind(element);
// // auto indicatorFunction = material_.getLocalIndicatorFunction();
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF);
// indicatorFunction.bind(element);
+ // Func2int indicatorFunction = *material_.getIndicatorFunction();
+
// Set of shape functions for a single element
const auto& localFiniteElement= localView.tree().child(0).finiteElement();
@@ -520,9 +585,13 @@ public:
defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
- double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
+ double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
+ // double energyDensity= scalarProduct(material_.ElasticityTensorPhase((-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));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
+ // double energyDensity= scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),defGradientV);
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal((-1.0)*forceTerm(quadPos),quadPos),defGradientV);
@@ -539,12 +608,16 @@ public:
// "f*m"-part
for (size_t m=0; m<3; m++)
{
- double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
+ double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
+ // double energyDensityfG= scalarProduct(material_.ElasticityTensorPhase((-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]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),basisContainer[m]);
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),basisContainer[m]);
+ // double energyDensityfG = scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),basisContainer[m]);
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 1fa1abe4f1b81d724f97bcacd145563ce8fcd6d6..90d09e54ca8a22915b8b012d644b0c408baf6a42 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -9,6 +9,8 @@
#include <dune/common/parametertree.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+
#include <dune/fufem/dunepython.hh>
@@ -24,10 +26,10 @@ using std::shared_ptr;
using std::make_shared;
-// constexpr unsigned int str2int(const std::string* str, int h = 0)
-// {
-// return !str[h] ? 5381 : (str2int(str, h+1) * 33) ^ str[h];
-// }
+
+
+
+
template<class Domain>
int indicatorFunction_material_1(const Domain& x)
@@ -61,8 +63,6 @@ using std::make_shared;
{
return 0;
}
-
-
MatrixRT material_homogeneous(const MatrixRT& G, const int& phase)
{
const FieldVector<double,1> mu = {80.0};
@@ -93,7 +93,19 @@ using std::make_shared;
-
+// template<class Domain>
+// MatrixRT MAT(const MatrixRT& G, const Domain& x) //unfortunately not possible as a GridViewFunction?
+// {
+// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
+// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
+// double theta=0.25;
+// if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
+// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
+// else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
+// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
+// else
+// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
+// }
@@ -105,10 +117,9 @@ public:
static const int dimworld = 3; //GridView::dimensionworld;
static const int dim = 3; //const int dim = Domain::dimension;
-
- // using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- // using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using LocalDomain = typename GridView::template Codim<0>::Geometry::LocalCoordinate;
+ using Element = typename GridViewEntitySet<GridView, 0>::Element;
using ScalarRT = FieldVector< double, 1>;
using VectorRT = FieldVector< double, dimworld>;
using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
@@ -119,7 +130,6 @@ public:
using Func2TensorPhase = std::function< MatrixRT(const MatrixRT& ,const int&) >;
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
-
protected:
const GridView& gridView_;
@@ -132,7 +142,6 @@ protected:
MatrixFunc L3_;
-
// const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
// const FuncScalar mu_;
@@ -153,78 +162,43 @@ protected:
// FuncScalar indicatorFunction_;
- // GridViewFunction<double(const Domain&), GridView> indicatorFunction_;
- GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
Func2int indicatorFunction_;
- // static const auto indicatorFunction_;
-
-// VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
-// VectorCT phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
-// FieldVector<double,3> m_1_, m_2_, m_3_; // Corrector m_i coefficient vectors
-
-// MatrixRT M1_, M2_, M3_; // (assembled) corrector matrices M_i
-// const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
-// const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
-
- // ---- Basis for R_sym^{2x2}
- MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
- MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G3_ {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
- std::array<MatrixRT,3 > MatrixBasisContainer_ = {G1_, G2_, G3_};
-
- Func2Tensor x3G_1_ = [] (const Domain& x) {
- return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
- };
-
- Func2Tensor x3G_2_ = [] (const Domain& x) {
- return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
- };
-
- Func2Tensor x3G_3_ = [] (const Domain& x) {
- return MatrixRT{{0.0, (1.0/sqrt(2.0))*x[2], 0.0}, {(1.0/sqrt(2.0))*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
-
- const std::array<Func2Tensor, 3> x3MatrixBasisContainer_ = {x3G_1_, x3G_2_, x3G_3_};
+ GridViewFunction<int(const Domain&), GridView> indicatorFunctionGVF_;
+ LocalFunction<int(const Domain&), typename GridViewEntitySet<GridView, 0>::Element > localIndicatorFunction_;
-
public:
///////////////////////////////
// constructor
///////////////////////////////
prestrainedMaterial(const GridView gridView,
- const ParameterTree& parameterSet) // string: "name of material"? // mu_(mu), muValues? müsste Anzahl Phasen bereits kennen..
+ const ParameterTree& parameterSet)
: gridView_(gridView),
parameterSet_(parameterSet)
{
std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
-
-
std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
// Python::Module module = Python::import(materialFunctionName_);
-
// elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
// elasticityTensor_ = setupElasticityTensor(materialFunctionName_);
setupElasticityTensor(materialFunctionName_);
- // elasticityTensor_ = material_1;
- // elasticityTensor_ = setupElasticityTensor();
- // module.get("Phases").toC<int>(Phases_);
+
+ indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
+ localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
// 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_);
// L1_ = Python::make_function<MatrixRT>(module.get("L1"));
// L2_ = Python::make_function<MatrixRT>(module.get("L2"));
// L3_ = Python::make_function<MatrixRT>(module.get("L3"));
-
-
- // Func2TensorParam elasticityTensor_ = Python::make_function<double>(module.get("L"));
- // Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
- // bool isotropic_ = true; // read from module File TODO
+ // bool isotropic_ = true; // read from module File
+ // auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
}
/////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -249,66 +223,6 @@ public:
// }
- // static int indicatorFunction_material_1(const Domain& x)
- // {
- // // std::cout << "x[0] , x[1] , x[2]" << x[0] << " , " << x[1] << ", "<< x[2] << std::endl;
- // double theta=0.25;
- // if ((x[0] < (-1/2+theta)) and (x[2]<(-1/2+theta)))
- // return 1; //#Phase1
- // else if ((x[1]< (-1/2+theta)) and (x[2]>(1/2-theta)))
- // return 2; //#Phase2
- // else
- // return 0; //#Phase3
- // }
-
- static MatrixRT material_1(const MatrixRT& G, const int& phase)
- {
- const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
- const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
- if (phase == 1)
- return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
- if (phase == 2)
- return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
- else
- return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
- }
-
-
-// static int indicatorFunction_material_1(const Domain& x)
-// {
-// double theta=0.25;
-// if (x[0] < (-1/2+theta) && x[2]<(-1/2+theta))
-// return 1; //#Phase1
-// else if (x[1]< (-1/2+theta) && x[2]>(1/2-theta))
-// return 2; //#Phase2
-// else
-// return 0; //#Phase3
-// }
-// static MatrixRT material_1(const MatrixRT& G, const int& phase)
-// {
-// FieldVector<double,3> mu = {80.0, 80.0, 60.0};
-// FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
-// if (phase == 1)
-// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
-// if (phase == 2)
-// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
-// else
-// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
-// }
-
-
-//---function that determines elasticity Tensor
- // auto setupElasticityTensor(const std::string name) // cant use materialFunctionName_ here!?
- // {
- // if(name == "material")
- // {
- // return material_1;
- // }
- // else
- // DUNE_THROW(Exception, "There exists no material in materialDefinitions.hh with this name ");
- // }
-
-
//---function that determines elasticity Tensor
void setupElasticityTensor(const std::string name) // cant use materialFunctionName_ here!?
@@ -318,7 +232,7 @@ public:
{
// indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// indicatorFunction_ = indicatorFunction_material_1;
- indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1<Domain>, gridView_);
+
indicatorFunction_ = indicatorFunction_material_1<Domain>;
elasticityTensor_ = material_1;
}
@@ -326,7 +240,7 @@ public:
{
// indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// indicatorFunction_ = indicatorFunction_material_1;
- indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_homogeneous<Domain>, gridView_);
+ // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_homogeneous<Domain>, gridView_);
indicatorFunction_ = indicatorFunction_material_homogeneous<Domain>;
elasticityTensor_ = material_homogeneous;
}
@@ -334,7 +248,7 @@ public:
{
// indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// indicatorFunction_ = indicatorFunction_material_1;
- indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_2<Domain>, gridView_);
+ // indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_material_2<Domain>, gridView_);
indicatorFunction_ = indicatorFunction_material_2<Domain>;
elasticityTensor_ = material_2;
}
@@ -343,10 +257,30 @@ public:
}
+ // MatrixRT localElasticityTensor(const MatrixRT& G, const Domain& x, Element& element) const //local Version ... muss binding öfters durchführen
+ // {
+ // localIndicatorFunction_.bind(*element);
+ // return elasticityTensor_(G,localIndicatorFunction_(x));
+ // }
+
+
+
+ //--- apply elasticityTensor_ to input Matrix G at position x
+ MatrixRT ElasticityTensorPhase(const MatrixRT& G, const int& phase) const //global Version (takes global coordinates)
+ {
+ return elasticityTensor_(G,phase);
+
+ }
+ MatrixRT ElasticityTensorGlobal(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
+ {
+ return MAT(G,x);
+
+ }
//--- apply elasticityTensor_ to input Matrix G at position x
- MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
+ MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
{
// Python::Module module = Python::import("material");
+
// auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
// return elasticityTensor_(G,indicatorFunctionTest(x));
// auto IFunction = localFunction(indicatorFunction_);
@@ -354,6 +288,7 @@ public:
// auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
// auto tmpLocal = LocalF
// return elasticityTensor_(G,IFunction(x));
+ // return elasticityTensor_(G,localIndicatorFunction_(x));
return elasticityTensor_(G,indicatorFunction_(x));
// return elasticityTensor_(G,indicatorFunctionGVF_(x));
// int phase = indicatorFunction_(x);
@@ -373,10 +308,6 @@ public:
//////////////////////////////////////////////////////////////////////////////
-
-
-
-
// MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
// {
// //--- apply elasticityTensor_ to input Matrix G at position x
@@ -410,74 +341,27 @@ public:
};
-
-
///////////////////////////////
// getter
///////////////////////////////
ParameterTree getParameterSet() const {return parameterSet_;}
-
-
// Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
// Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
-
-
- // auto getIndicatorFunction() const {return indicatorFunction_;}
- auto getLocalIndicatorFunction() const {return localFunction(indicatorFunctionGVF_);} //get as localFunction
+ auto getLocalIndicatorFunction() const {return localIndicatorFunction_;} //get as localFunction
auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;} //get as GridViewFunction
auto getIndicatorFunction() const {return indicatorFunction_;}
// shared_ptr<Func2int> getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
// auto getIndicatorFunction(){return make_shared<Func2int>(indicatorFunction_);}
-
-
-
-
// shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}
-
//TODO getLameParameters() .. Throw Exception if isotropic_ = false!
- // shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
- // shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
- // shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
- // shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
-
- // shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
- // shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
-
-
- // --- Get Correctors
- // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
- // auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
- // auto getMcontainer(){return mContainer;}
- // shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}
-
-
- // // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
- // auto getMatrixBasiscontainer(){return make_shared<std::array<MatrixRT,3 >>(MatrixBasisContainer_);}
- // // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
- // auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
-
-
-
-
-
- // shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
- // shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
- // shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
- // shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
-
-
-
-
-
-
}; // end class
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 3dcbfe0500fa7148cb14756e22417478f51d8a26..46f1ffae0eba37b1b1cac55ac4c8f108e3382992 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -25,7 +25,7 @@ print_debug = true #(default=false)
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
#----------------------------------------------------
-numLevels= 2 3
+numLevels= 2 4
#numLevels = 0 0 # computes all levels from first to second entry
#numLevels = 2 2 # computes all levels from first to second entry
#numLevels = 1 3 # computes all levels from first to second entry
@@ -47,7 +47,7 @@ gamma=1.0
#############################################
#-- stiffness-ratio (ratio between material parameters mu1 & mu2 .... beta = 1.0 corresponds to homogeneous case)
-beta=3.0
+beta=2.0
$--- strength of prestrain
rho1=1.0
@@ -57,8 +57,8 @@ alpha=8.0
#materialFunction = "homogeneous"
-#materialFunction = "material"
-materialFunction = "material_2"
+materialFunction = "material_1"
+#materialFunction = "material_2"
#--- Lame-Parameters
mu1=80.0
@@ -90,10 +90,10 @@ theta=0.25
#--- choose composite-Implementation:
#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
geometryFunctionPath = /home/klaus/Desktop/Dune-Testing/dune-microstructure/geometries
-#material_prestrain_imp= "material_neukamm" #(Python-indicator-function with same name determines material phases)
+material_prestrain_imp= "material_neukamm" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_2" #(Python-indicator-function with same name determines material phases)
#material_prestrain_imp= "two_phase_material_3" #(Python-indicator-function with same name determines material phases)
-material_prestrain_imp= "parametrized_Laminate"
+#material_prestrain_imp= "parametrized_Laminate"
#material_prestrain_imp= "homogeneous"
#material_prestrain_imp= "analytical_Example"
#material_prestrain_imp="isotropic_bilayer"
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 537c9c4443ac8180445d386eb999d563d21353b9..fdbf71640e52166182508a0f886fa11116449dda 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -42,7 +42,7 @@
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrix.hh>
#include <dune/microstructure/prestrain_material_geometry.hh>
#include <dune/microstructure/matrix_operations.hh>