diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 762125683373d21057d427b7c73397f7168da99e..7aef0e24ad7ceb45a3f0ddd3508627afc1b2e944 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -2,7 +2,7 @@
#define DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#include <dune/common/parametertree.hh>
-#include <dune/common/float_cmp.hh>
+// #include <dune/common/float_cmp.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
@@ -13,7 +13,6 @@
#include <dune/solvers/solvers/umfpacksolver.hh>
using namespace Dune;
-// using namespace Functions;
using namespace MatrixOperations;
using std::shared_ptr;
@@ -30,16 +29,13 @@ public:
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-
using ScalarRT = FieldVector< double, 1>;
using VectorRT = FieldVector< double, dimworld>;
using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
-
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> >;
using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
@@ -56,8 +52,8 @@ protected:
fstream& log_; // Output-log
const ParameterTree& parameterSet_;
- const FuncScalar mu_;
- const FuncScalar lambda_;
+ // const FuncScalar mu_;
+ // const FuncScalar lambda_;
double gamma_;
MatrixCT stiffnessMatrix_;
@@ -100,15 +96,11 @@ public:
///////////////////////////////
CorrectorComputer( const Basis& basis,
const Material& material,
- const FuncScalar& mu,
- const FuncScalar& lambda,
double gamma,
std::fstream& log,
const ParameterTree& parameterSet)
: basis_(basis),
material_(material),
- mu_(mu),
- lambda_(lambda),
gamma_(gamma),
log_(log),
parameterSet_(parameterSet),
@@ -117,12 +109,7 @@ public:
assemble();
- // if (parameterSet.get<bool>("stiffnessmatrix_cellproblem_to_csv"))
- // csvSystemMatrix();
- // if (parameterSet.get<bool>("rhs_cellproblem_to_csv"))
- // csvRHSs();
- // if (parameterSet.get<bool>("rhs_cellproblem_to_vtk"))
- // vtkLoads();
+ //write VTK call here...
}
@@ -157,8 +144,8 @@ public:
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_);}
+ // shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
+ // shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
@@ -174,9 +161,6 @@ public:
// 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_);}
@@ -249,12 +233,15 @@ public:
std::cout << "finished occupation pattern\n";
}
+ // template<class localFunction1, class localFunction2>
+ // void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
+ // ElementMatrixCT& elementMatrix,
+ // const localFunction1& mu,
+ // const localFunction2& lambda
+ // )
- template<class localFunction1, class localFunction2>
void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
- ElementMatrixCT& elementMatrix,
- const localFunction1& mu,
- const localFunction2& lambda
+ ElementMatrixCT& elementMatrix
)
{
// Local StiffnessMatrix of the form:
@@ -485,10 +472,15 @@ public:
// Compute the source term for a single element
// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
- template<class LocalFunction1, class LocalFunction2, class Vector, class LocalForce>
+ // template<class LocalFunction1, class LocalFunction2, class Vector, class LocalForce>
+ // void computeElementLoadVector( const typename Basis::LocalView& localView,
+ // LocalFunction1& mu,
+ // LocalFunction2& lambda,
+ // Vector& elementRhs,
+ // const LocalForce& forceTerm
+ // )
+ template<class Vector, class LocalForce>
void computeElementLoadVector( const typename Basis::LocalView& localView,
- LocalFunction1& mu,
- LocalFunction2& lambda,
Vector& elementRhs,
const LocalForce& forceTerm
)
@@ -641,21 +633,23 @@ public:
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
- auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- auto muLocal = localFunction(muGridF);
- auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- auto lambdaLocal = localFunction(lambdaGridF);
+ // auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ // auto muLocal = localFunction(muGridF);
+ // auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ // auto lambdaLocal = localFunction(lambdaGridF);
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
- muLocal.bind(element);
- lambdaLocal.bind(element);
+ // muLocal.bind(element);
+ // lambdaLocal.bind(element);
const int localPhiOffset = localView.size();
// Dune::Timer Time;
//std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
Matrix<FieldMatrix<double,1,1> > elementMatrix;
- computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
+ // computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
+ computeElementStiffnessMatrix(localView, elementMatrix);
+
// std::cout << "local assembly time:" << Time.elapsed() << std::endl;
@@ -712,18 +706,18 @@ public:
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis_.gridView());
auto loadFunctional = localFunction(loadGVF);
- auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- auto muLocal = localFunction(muGridF);
- auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- auto lambdaLocal = localFunction(lambdaGridF);
+ // auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ // auto muLocal = localFunction(muGridF);
+ // auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ // auto lambdaLocal = localFunction(lambdaGridF);
// int counter = 1;
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
- muLocal.bind(element);
- lambdaLocal.bind(element);
+ // muLocal.bind(element);
+ // lambdaLocal.bind(element);
loadFunctional.bind(element);
const int localPhiOffset = localView.size();
@@ -732,7 +726,8 @@ public:
VectorCT elementRhs;
// std::cout << "----------------------------------Element : " << counter << std::endl; //TEST
// counter++;
- computeElementLoadVector(localView, muLocal, lambdaLocal, elementRhs, loadFunctional);
+ // computeElementLoadVector(localView, muLocal, lambdaLocal, elementRhs, loadFunctional);
+ computeElementLoadVector(localView, elementRhs, loadFunctional);
// computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
// printvector(std::cout, elementRhs, "elementRhs", "--");
// printvector(std::cout, elementRhs, "elementRhs", "--");
@@ -1304,10 +1299,10 @@ public:
- auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- auto mu = localFunction(muGridF);
- auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- auto lambda= localFunction(lambdaGridF);
+ // auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ // auto mu = localFunction(muGridF);
+ // auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ // auto lambda= localFunction(lambdaGridF);
for(int a=0; a<3; a++)
for(int b=0; b<3; b++)
@@ -1339,8 +1334,8 @@ public:
matrixFieldG.bind(e);
DerPhi1.bind(e);
DerPhi2.bind(e);
- mu.bind(e);
- lambda.bind(e);
+ // mu.bind(e);
+ // lambda.bind(e);
localIndicatorFunction.bind(e);
double elementEnergy = 0.0;
@@ -1436,8 +1431,6 @@ public:
std::cout << "--- Symmetry test passed ---" << std::endl;
}
-
-
}; // end class
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index c7f6add86a71ea4cfbab5ea34a3c52cfe7be85ac..7dcb8d2bbb76bb4aa2735d3c7aebceea5f6721f3 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -1,12 +1,9 @@
#ifndef DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
-#include <filesystem>
-
-
+// #include <filesystem>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/CorrectorComputer.hh>
-
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
#include <dune/istl/io.hh>
#include <dune/istl/matrix.hh>
@@ -29,8 +26,6 @@ class EffectiveQuantitiesComputer {
public:
static const int dimworld = 3;
static const int dim = Basis::GridView::dimension;
-
-
using Domain = typename CorrectorComputer<Basis,Material>::Domain;
using VectorRT = typename CorrectorComputer<Basis,Material>::VectorRT;
using MatrixRT = typename CorrectorComputer<Basis,Material>::MatrixRT;
@@ -42,93 +37,41 @@ public:
using MatrixPhase = std::function< MatrixRT(const int&) >;
protected:
-
CorrectorComputer<Basis,Material>& correctorComputer_;
// Func2Tensor prestrain_;
MatrixPhase prestrain_;
const Material& material_;
-
// Func2int indicatorFunction_;
public:
- VectorCT B_load_TorusCV_; //<B, Chi>_L2
- // FieldMatrix<double, dim, dim> Q_; //effective moduli <LF_i, F_j>_L2
- // FieldVector<double, dim> Bhat_; //effective loads induced by prestrain <LF_i, B>_L2
- // FieldVector<double, dim> Beff_; //effective strains Mb = ak
- MatrixRT Q_; //effective moduli <LF_i, F_j>_L2
- VectorRT Bhat_; //effective loads induced by prestrain <LF_i, B>_L2
- VectorRT Beff_; //effective strains Mb = ak
+
+ MatrixRT Q_; //effective moduli
+ VectorRT Bhat_;
+ VectorRT Beff_;
// corrector parts
- VectorCT phi_E_TorusCV_; //phi_i * (a,K)_i
- VectorCT phi_perp_TorusCV_;
- VectorCT phi_TorusCV_;
- VectorCT phi_1_; //phi_i * (a,K)_i
- VectorCT phi_2_;
- VectorCT phi_3_;
+ // VectorCT phi_1_;
+ // VectorCT phi_2_;
+ // VectorCT phi_3_;
- // is this really interesting???
- // double phi_E_L2norm_;
- // double phi_E_H1seminorm_;
-
- // double phi_perp_L2norm_;
- // double phi_perp_H1seminorm_;
-
- // double phi_L2norm_;
- // double phi_H1seminorm_;
-
- // double Chi_E_L2norm_;
- // double Chi_perp_L2norm_;
- // double Chi_L2norm_;
-
-
- double B_energy_; // < B, B >_L B = F + Chi_perp + B_perp
- double F_energy_; // < F, F >_L
- double Chi_perp_energy_; // < Chi_perp, Chi_perp >_L
- double B_perp_energy_; // < B_perp, B_perp >_L
-
- //Chi(phi) is only implicit computed, can we store this?
///////////////////////////////
// constructor
///////////////////////////////
- // EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer, Func2Tensor prestrain)
- // : correctorComputer_(correctorComputer), prestrain_(prestrain)
- // EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
- // Func2Tensor prestrain,
- // const Material& material)
- // : correctorComputer_(correctorComputer),
- // prestrain_(prestrain),
- // material_(material)
EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
const Material& material)
: correctorComputer_(correctorComputer),
material_(material)
{
-
-
prestrain_ = material_.getPrestrainFunction();
- // computePrestressLoadCV();
- // computeEffectiveStrains();
- // Q_ = 0;
- // Q_ = {{0.0,0.0,0.0},{0.0,0.0,0.0},{0.0,0.0,0.0}};
- // compute_phi_E_TorusCV();
- // compute_phi_perp_TorusCV();
- // compute_phi_TorusCV();
-
- // computeCorrectorNorms();
- // computeChiNorms();
- // computeEnergiesPrestainParts();
-
- // writeInLogfile();
}
///////////////////////////////
- // getter
+ // Getter
///////////////////////////////
CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}
@@ -152,8 +95,8 @@ public:
auto MContainer = correctorComputer_.getMcontainer();
auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
- auto mu_ = *correctorComputer_.getMu();
- auto lambda_ = *correctorComputer_.getLambda();
+ // auto mu_ = *correctorComputer_.getMu();
+ // auto lambda_ = *correctorComputer_.getLambda();
auto gamma = correctorComputer_.getGamma();
auto basis = *correctorComputer_.getBasis();
ParameterTree parameterSet = correctorComputer_.getParameterSet();
@@ -191,10 +134,10 @@ public:
auto matrixFieldG2GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[b], basis.gridView());
auto matrixFieldG2 = localFunction(matrixFieldG2GVF);
- auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
- auto mu = localFunction(muGridF);
- auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
- auto lambda= localFunction(lambdaGridF);
+ // auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
+ // auto mu = localFunction(muGridF);
+ // auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
+ // auto lambda= localFunction(lambdaGridF);
// using GridView = typename Basis::GridView;
@@ -205,8 +148,8 @@ public:
matrixFieldG2.bind(e);
localfun_a.bind(e);
// DerPhi2.bind(e);
- mu.bind(e);
- lambda.bind(e);
+ // mu.bind(e);
+ // lambda.bind(e);
// prestrainFunctional.bind(e);
localIndicatorFunction.bind(e);
@@ -303,9 +246,6 @@ public:
return;
}
-
-
-
template<class MatrixFunction>
double energySP(const MatrixFunction& matrixFieldFuncA,
const MatrixFunction& matrixFieldFuncB)
@@ -317,21 +257,24 @@ public:
auto basis = *correctorComputer_.getBasis();
auto localView = basis.localView();
+ auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+
auto matrixFieldAGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncA, basis.gridView());
auto matrixFieldA = localFunction(matrixFieldAGVF);
auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
auto matrixFieldB = localFunction(matrixFieldBGVF);
- auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
- auto mu = localFunction(muGridF);
- auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
- auto lambda= localFunction(lambdaGridF);
+ // auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
+ // auto mu = localFunction(muGridF);
+ // auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
+ // auto lambda= localFunction(lambdaGridF);
for (const auto& e : elements(basis.gridView()))
{
localView.bind(e);
matrixFieldA.bind(e);
matrixFieldB.bind(e);
- mu.bind(e);
- lambda.bind(e);
+ // mu.bind(e);
+ // lambda.bind(e);
+ localIndicatorFunction.bind(e);
double elementEnergy = 0.0;
@@ -344,7 +287,9 @@ public:
{
const auto& quadPos = quadPoint.position();
const double integrationElement = geometry.integrationElement(quadPos);
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), matrixFieldA(quadPos), matrixFieldB(quadPos));
+
+ double energyDensity= scalarProduct(material_.ElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), matrixFieldA(quadPos), matrixFieldB(quadPos));
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index afd2dca574b83bf62a5e9fcb59882a711e922fbb..9ce64e42022b1f616825ee17695e6a3c70a6fd8b 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -1,7 +1,6 @@
#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
-
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
@@ -24,31 +23,14 @@ using std::shared_ptr;
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();
-// }
-
-template <class GridView> // needed for GridViewFunctions
+//--------------------------------------------------------
+template <class GridView>
class prestrainedMaterial
{
public:
- static const int dimworld = 3; //GridView::dimensionworld;
- static const int dim = 3; //const int dim = Domain::dimension;
-
+ static const int dimworld = 3; //GridView::dimensionworld;
+ static const int dim = 3; //Domain::dimension;
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;
@@ -149,12 +131,11 @@ public:
}
- //--- apply elasticityTensor_ to input Matrix G at position x
+ //--- Apply elasticityTensor_ to input Matrix G at material phase
// input: Matrix G, material-phase
MatrixRT ElasticityTensor(const MatrixRT& G, const int& phase) const
{
return elasticityTensor_(G,phase);
-
}
int IndicatorFunction(const Domain& x) const
@@ -162,88 +143,6 @@ public:
return indicatorFunction_(x);
}
- // 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;
- // std::cout << "-1/2" << -1/2 << std::endl;
- // if ((x[1]< (-0.5+theta)) and (x[2]>(0.5-theta)))
- // {
- // return 2; //#Phase1
- // }
- // else if ((x[0] < (-0.5+theta)) and (x[2]<(-0.5+theta)))
- // {
- // return 1; //#Phase2
- // }
- // else
- // {
- // return 0; //#Phase3
- // }
- // }
-
-
-
- // 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));
- // }
-
-
-
-
- // 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 //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_);
- // // auto IFunction = this->getIndicatorFunction();
- // // 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);
- // // auto tmp = elasticityTensor_(G,phase);
- // // auto tmp = material_1(G,indicatorFunction_(x));
- // // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
- // // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
- // // return tmp;
- // // return material_1(G,indicatorFunction_(x));
- // }
-
-
-
-//////////////////////////////////////////////////////////////////////////////
- // MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
- // {
- // //--- apply elasticityTensor_ to input Matrix G at position x
- // return elasticityTensor_(G,x);
-
- // }
-
- // MatrixRT applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const
- // {
- // //--- apply elasticityTensor_ to input Matrix G at position x (local coordinates)
- // // MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
-
- // if (indicatorFunction_(x) == 1)
- // return L1_(G);
- // else if (indicatorFunction_(x) == 2)
- // return L2_(G);
- // else
- // return L3_(G);
-
- // }
///////////////////////////////
// Getter
@@ -258,13 +157,7 @@ public:
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_);}
-
- //getLameParameters() .. Throw Exception if isotropic_ = false!
-
+ //getLameParameters()? .. Throw Exception if isotropic_ = false!
///////////////////////////////
// VTK - Writer
@@ -380,6 +273,48 @@ public:
// // };
+ // 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 //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_);
+ // // auto IFunction = this->getIndicatorFunction();
+ // // 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);
+ // // auto tmp = elasticityTensor_(G,phase);
+ // // auto tmp = material_1(G,indicatorFunction_(x));
+ // // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
+ // // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
+ // }
+
+
+// template<class Domain>
+// static 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();
+// }
+
+
}; // end class
diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 46160e703402550e1d8a8bbb172f5ec601b61a61..7db3997ed4d640f52e737e9b06be62b938657987 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -34,7 +34,7 @@ numLevels= 2 2 # 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_2"
#materialFunction = "homogeneous"
@@ -43,12 +43,6 @@ gamma=1.0
-mu=80 80 60
-lambda=80 80 25
-rho1=1.0
-beta=3.0
-theta=0.25
-material_prestrain_imp= "material_neukamm"
#--- choose composite-Implementation:
#geometryFunctionPath = /home/stefan/DUNE/dune-microstructure/geometries
@@ -80,8 +74,8 @@ Solver_verbosity = 0 #(default = 2) degree of information for solver output
#############################################
# --- (Optional output) write Material / prestrain / Corrector functions to .vtk-Files:
-write_materialFunctions = true
-write_prestrainFunctions = true # VTK norm of B ,
+write_materialFunctions = true # VTK indicator function for material/prestrain definition
+#write_prestrainFunctions = true # VTK norm of B (currently not implemented)
# --- Write Correctos to VTK-File:
write_VTK = true
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index ad05b62ea2d77c6fe18c34994767d5aff60f2683..b639fb8b407fbb817e50789a37c66031829ba225 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -113,8 +113,6 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
-// // a function:
-// int half(int x, int y) {return x/2+y/2;}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -145,14 +143,14 @@ int main(int argc, char *argv[])
// parameterSet.report(log); // short Alternativ
//--- Get Path for Material/Geometry functions
- std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath");
+ // 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 << "sys.path.append('" << geometryFunctionPath << "')"
<< std::endl;
@@ -166,56 +164,6 @@ int main(int argc, char *argv[])
bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false);
-
-
-
- ///////////////////////////////////
- // Get Parameters/Data
- ///////////////////////////////////
- double gamma = parameterSet.get<double>("gamma",1.0); // ratio dimension reduction to homogenization
- double alpha = parameterSet.get<double>("alpha", 2.0);
- double theta = parameterSet.get<double>("theta",1.0/4.0);
- ///////////////////////////////////
- // Get Material Parameters
- ///////////////////////////////////
- std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
- log << "material_prestrain used: "<< imp << std::endl;
- double beta = parameterSet.get<double>("beta",2.0);
- double mu1 = parameterSet.get<double>("mu1",1.0);;
- double mu2 = beta*mu1;
- double lambda1 = parameterSet.get<double>("lambda1",0.0);;
- double lambda2 = beta*lambda1;
-
-
- if(imp == "material_neukamm")
- {
- std::cout <<"mu: " << parameterSet.get<std::array<double,3>>("mu", {1.0,3.0,2.0}) << std::endl;
- std::cout <<"lambda: " << parameterSet.get<std::array<double,3>>("lambda", {1.0,3.0,2.0}) << std::endl;
- }
- else
- {
- std::cout <<"mu: " << parameterSet.get<double>("mu1",1.0) << std::endl;
- std::cout <<"lambda: " << parameterSet.get<double>("lambda1",0.0) << std::endl;
- }
-
- ///////////////////////////////////
- // Get Prestrain/Parameters
- ///////////////////////////////////
- auto prestrainImp = PrestrainImp<dim>();
- auto B_Term = prestrainImp.getPrestrain(parameterSet);
-
-
-
- log << "----- Input Parameters -----: " << std::endl;
- log << "alpha: " << alpha << std::endl;
- log << "gamma: " << gamma << std::endl;
- log << "theta: " << theta << std::endl;
- log << "beta: " << beta << std::endl;
- log << "material parameters: " << std::endl;
- log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
- log << "lambda1: " << lambda1 <<"\nlambda2: " << lambda2 << std::endl;
- log << "----------------------------: " << std::endl;
-
///////////////////////////////////
// Generate the grid
///////////////////////////////////
@@ -228,15 +176,16 @@ int main(int argc, char *argv[])
///////////////////////////////////
- // Create Data Storage
+ // Create Date 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_Quantities;
- // for(const size_t &level : numLevels) // explixite Angabe der levels.. {2,4}
- for(size_t level = numLevels[0] ; level <= numLevels[1]; level++) // levels von bis.. [2,4]
+
+ //--- GridLevel-Loop:
+ for(size_t level = numLevels[0] ; level <= numLevels[1]; level++)
{
std::cout << " ----------------------------------" << std::endl;
std::cout << "GridLevel: " << level << std::endl;
@@ -254,27 +203,6 @@ int main(int argc, char *argv[])
const GridView gridView_CE = grid_CE.leafGridView();
if(print_debug)
std::cout << "Host grid has " << gridView_CE.size(dim) << " vertices." << std::endl;
-
- // //not needed
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
- using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
- using Func2Tensor = std::function< MatrixRT(const Domain&) >;
- // using Func2Tensor = std::function< MatrixRT(const Domain&) >;
- // using VectorCT = BlockVector<FieldVector<double,1> >;
- // using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
-
- ///////////////////////////////////
- // Create Lambda-Functions for material Parameters depending on microstructure
- ///////////////////////////////////
- auto materialImp = IsotropicMaterialImp<dim>();
- auto muTerm = materialImp.getMu(parameterSet);
- auto lambdaTerm = materialImp.getLambda(parameterSet);
-
- auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
- auto muLocal = localFunction(muGridF);
- auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
- auto lambdaLocal = localFunction(lambdaGridF);
-
//--- Choose a finite element space for Cell Problem
using namespace Functions::BasisFactory;
@@ -300,293 +228,34 @@ int main(int argc, char *argv[])
std::cout << "power<periodic> basis has " << Basis_CE.dimension() << " degrees of freedom" << std::endl;
-
-
-
- //TEST
- //Read from Parset...
- // int Phases = parameterSet.get<int>("Phases", 3);
-
-
- // std::string materialFunctionName = parameterSet.get<std::string>("materialFunction", "material");
- // Python::Module module = Python::import(materialFunctionName);
-
-
-
- // auto indicatorFunction = Python::make_function<double>(module.get("f"));
- // Func2Tensor indicatorFunction = Python::make_function<double>(module.get("f"));
- // auto materialFunction_ = Python::make_function<double>(module.get("f"));
- // auto materialFunction_ = Python::make_function<double>(module.get("f"));
- // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
-
- // int Phases;
- // module.get("Phases").toC<int>(Phases);
- // std::cout << "Number of Phases used:" << Phases << std::endl;
-
-
- // std::cout << typeid(mu_).name() << '\n';
-
-
- //---- Get mu/lambda values (for isotropic material) from Material-file
- // FieldVector<double,3> mu_(0);
- // module.get("mu_").toC<FieldVector<double,3>>(mu_);
- // printvector(std::cout, mu_ , "mu_", "--");
- // FieldVector<double,3> lambda_(0);
- // module.get("lambda_").toC<FieldVector<double,3>>(lambda_);
- // printvector(std::cout, lambda_ , "lambda_", "--");
-
-
- //////////////////////////////////////////////////////////////////////////////////////////////////////////
- // TESTING PRESTRAINEDMATERIAL.HH:
- using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
-
+ ///////////////////////////////////
+ // Create prestrained material object
+ ///////////////////////////////////
auto material_ = prestrainedMaterial(gridView_CE,parameterSet);
- // Func2Tensor elasticityTensor = material_.getElasticityTensor();
- // auto elasticityTensor = material_.getElasticityTensor();
- // Func2TensorParam elasticityTensor_ = *material_.getElasticityTensor();
- auto elasticityTensor_ = material_.getElasticityTensor();
-
- // Func2TensorParam TestTensor = Python::make_function<MatrixRT>(module.get("H"));
-
- // std::cout << "decltype(elasticityTensor_) " << decltype(elasticityTensor_) << std::endl;
-
-
- // std::cout << "typeid(elasticityTensor).name() :" << typeid(elasticityTensor_).name() << '\n';
- // std::cout << "typeid(TestTensor).name() :" << typeid(TestTensor).name() << '\n';
-
- using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
- // std::cout << "Import NOW:" << std::endl;
- // MatrixFunc symTest = Python::make_function<MatrixRT>(module.get("sym"));
- // auto indicatorFunction = Python::make_function<double>(module.get("indicatorFunction"));
-
- //localize..
- // auto indicatorFunctionGVF = Dune::Functions::makeGridViewFunction(indicatorFunction , Basis_CE.gridView());
- // GridViewFunction<double(const Domain&), GridView> indicatorFunctionGVF = Dune::Functions::makeGridViewFunction(indicatorFunction , Basis_CE.gridView());
- // auto localindicatorFunction = localFunction(indicatorFunctionGVF);
-
-
- // auto indicatorFunctionGVF = material_.getIndicatorFunction();
- // auto localindicatorFunction = localFunction(indicatorFunctionGVF);
-
- auto localindicatorFunction = material_.getIndicatorFunction();
-
- // GridView::Element
- // auto localindicatorFunction = localFunction(indicatorFunction);
-
- // std::cout << "typeid(localindicatorFunction).name() :" << typeid(localindicatorFunction).name() << '\n';
-
- // using MatrixDomainFunc = std::function< MatrixRT(const MatrixRT&,const Domain&)>;
- // // MatrixFunc elasticityTensor = Python::make_function<MatrixRT>(module.get("L"));
-
- // auto elasticityTensorGVF = Dune::Functions::makeGridViewFunction(elasticityTensor , Basis_CE.gridView());
- // auto localElasticityTensor = localFunction(elasticityTensorGVF);
-
- // Func2Tensor forceTerm = [] (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?
- // };
-
- // auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, Basis_CE.gridView());
- // auto loadFunctional = localFunction(loadGVF);
-
- MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
- // auto xu = symTest(G1_);
- // std::cout << "TEST NOW:" << std::endl;
- // printmatrix(std::cout, symTest(G1_), "symTest(G1_)", "--");
-
- // auto TestTensorGVF = Dune::Functions::makeGridViewFunction(TestTensor , Basis_CE.gridView());
- // auto localTestTensor = localFunction(TestTensorGVF );
-
-
- // printmatrix(std::cout, elasticityTensor(G1_), "elasticityTensor(G1_)", "--");
- // auto temp = elasticityTensor(G1_);
-
-
-
- // for (const auto& element : elements(Basis_CE.gridView()))
- // {
- // localindicatorFunction.bind(element);
- // int orderQR = 2;
- // const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
- // for (const auto& quadPoint : quad)
- // {
- // const auto& quadPos = quadPoint.position();
-
- // // std::cout << "localindicatorFunction(quadPos): " << localindicatorFunction(quadPos) << std::endl;
-
-
-
- // // std::cout << "quadPos : " << quadPos << std::endl;
- // // auto temp = TestTensor(G1_, element.geometry().global(quadPos));
- // // auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
- // // std::cout << "material_.applyElasticityTensor:" << std::endl;
- // // auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos));
- // // auto tmp3 = material_.applyElasticityTensor(G1_, quadPos);
- // // printmatrix(std::cout, tmp3, "tmp3", "--");
- // }
- // }
-
-
-
- // for (auto&& vertex : vertices(gridView_CE))
- // {
- // std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
- // auto tmp = vertex.geometry().corner(0);
- // auto temp = elasticityTensor(tmp);
- // // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
- // // printmatrix(std::cout, localElasticityTensor(G1_,tmp), "localElasticityTensor(vertex.geometry().corner(0))", "--");
- // }
-
-
-
- // std::function<int(int,int)> fn1 = half;
- // std::cout << "fn1(60,20): " << fn1(60,20) << '\n';
-
-
- // std::cout << typeid(elasticityTensorGVF).name() << '\n';
- // std::cout << typeid(localElasticityTensor).name() << '\n';
-
-
- // ParameterTree parameterSet_2;
- // ParameterTreeParser::readINITree(geometryFunctionPath + "/"+ materialFunctionName + ".py", parameterSet_2);
-
- // auto lu = parameterSet_2.get<FieldVector<double,3>>("lu", {1.0,3.0,2.0});
- // std::cout <<"lu[1]: " << lu[1]<< std::endl;
- // std::cout <<"lu: " << parameterSet_2.get<std::array<double,3>>("lu", {1.0,3.0,2.0}) << std::endl;
-
- // auto mU_ = module.evaluate(parameterSet_2.get<std::string>("lu", "[1,2,3]"));
- // std::cout << "typeid(mU_).name()" << typeid(mU_.operator()()).name() << '\n';
-
-
- // for (auto&& vertex : vertices(gridView_CE))
- // {
- // std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
- // // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
- // printvector(std::cout, materialFunction_(vertex.geometry().corner(0)), "materialFunction_(vertex.geometry().corner(0))", "--");
- // }
- // std::cout << "materialFunction_({0.0,0.0,0.0})", materialFunction_({0.0,0.0,0.0}) << std::endl;
-
-
-
- // --------------------------------------------------------------
-
- //TODO// Phasen anhand von Mu bestimmen?
- //TODO: DUNE_THROW(Exception, "Inconsistent choice of interpolation method"); if number of Phases != mu/lambda parameters
-
- //FÜR L GARNICHT NÖTIG DENN RÜCKGABETYPE IS IMMER MATRIXRT!?!:
- // BEi materialfunction (isotopic) reicht auch FieldVector<double,2> für lambda/mu
-
- // switch (Phases)
- // {
- // case 1: //homogeneous material
- // {
- // std::cout << "Phase - 1" << std::endl;
- // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
- // break;
- // }
- // case 2:
- // {
- // std::cout << "Phase - 1" << std::endl;
- // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
- // break;
- // }
- // case 3:
- // {
- // std::cout << "Phase - 3" << std::endl;
- // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
- // break;
- // }
- // }
-
-
- // switch (Phases)
- // {
- // case 1: //homogeneous material
- // {
- // std::cout << "Phases - 1" << std::endl;
- // std::array<double,1> mu_ = parameterSet.get<std::array<double,1>>("mu", {1.0});
- // Python::Module module = Python::import(materialFunction);
- // auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
- // auto muTerm = [mu_] (const Domain& x) {return mu_;};
- // break;
- // }
- // case 2:
- // {
- // std::cout << "Phases - 2" << std::endl;
- // std::array<double,2> mu_ = parameterSet.get<std::array<double,2>>("mu", {1.0,3.0});
- // Python::Module module = Python::import(materialFunction);
- // auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
- // auto muTerm = [mu_,indicatorFunction] (const Domain& x)
- // {
- // if (indicatorFunction(x) == 1)
- // return mu_[0];
- // else
- // return mu_[1];
- // };
- // break;
- // }
- // case 3:
- // {
- // std::cout << "Phases - 3" << std::endl;
- // std::array<double,3> mu_ = parameterSet.get<std::array<double,3>>("mu", {1.0,3.0, 5.0});
- // Python::Module module = Python::import(materialFunction);
- // auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
- // auto muTerm = [mu_,indicatorFunction] (const Domain& x)
- // {
- // if (indicatorFunction(x) == 1)
- // return mu_[0];
- // else if (indicatorFunction(x) == 2)
- // return mu_[1];
- // else
- // return mu_[2];
- // };
- // break;
- // }
- // }
-
-
-
-
- //TEST
-// std::cout << "Test crossSectionDirectionScaling:" << std::endl;
-/*
- MatrixRT T {{1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}};
- printmatrix(std::cout, T, "Matrix T", "--");
-
- auto ST = crossSectionDirectionScaling((1.0/5.0),T);
- printmatrix(std::cout, ST, "scaled Matrix T", "--");*/
-
- //TEST
-// auto QuadraticForm = [] (const double mu, const double lambda, const MatrixRT& M) {
-//
-// return lambda*std::pow(trace(M),2) + 2*mu*pow(norm(sym(M)),2);
-// };
-
-
+ // --- get scale ratio
+ double gamma = parameterSet.get<double>("gamma",1.0);
//------------------------------------------------------------------------------------------------
- //--- compute Correctors
+ //--- 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_, muTerm, lambdaTerm, gamma, log, parameterSet);
+ auto correctorComputer = CorrectorComputer(Basis_CE, material_, gamma, log, parameterSet);
correctorComputer.solve();
- //--- check Correctors (options):
+ //--- Check Correctors (options):
if(parameterSet.get<bool>("write_L2Error", false))
correctorComputer.computeNorms();
if(parameterSet.get<bool>("write_VTK", false))
correctorComputer.writeCorrectorsVTK(level);
- //--- additional Test: check orthogonality (75) from paper:
+ //--- Additional Test: check orthogonality (75) from paper:
if(parameterSet.get<bool>("write_checkOrthogonality", false))
correctorComputer.check_Orthogonality();
//--- Check symmetry of stiffness matrix
if(print_debug)
correctorComputer.checkSymmetry();
- // auto B_Term = material_.getPrestrainFunction();
-
-
- //--- compute effective quantities
- // auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term,material_);
+ //--- Compute effective quantities
auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material_);
effectiveQuantitiesComputer.computeEffectiveQuantities();
@@ -596,9 +265,7 @@ int main(int argc, char *argv[])
material_.writeVTKMaterialFunctions(nElements,level);
}
-
-
- //--- Test:: Compute Qeff without using the orthogonality (75)...
+ //--- TEST:: Compute Qeff without using the orthogonality (75)...
// only really makes a difference whenever the orthogonality is not satisfied!
// std::cout << "----------computeFullQ-----------"<< std::endl; //TEST
// effectiveQuantitiesComputer.computeFullQ();
@@ -620,15 +287,7 @@ int main(int argc, char *argv[])
std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Qeff[0][1] << std::endl;
// -------------------------------------------
-
- //TEST
- // 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?
- // };
-
- // double energy = effectiveQuantitiesComputer.energySP(x3G_1,x3G_1);
- // std::cout << "energy:" << energy << std::endl;
-
+ // --- Fill output-Table:
Storage_Quantities.push_back(Qeff[0][0] );
Storage_Quantities.push_back(Qeff[1][1] );
Storage_Quantities.push_back(Qeff[2][2] );
@@ -651,18 +310,12 @@ int main(int argc, char *argv[])
log << "mu_gamma=" << Qeff[2][2] << std::endl; // added for Python-Script
-
-
-
-
levelCounter++;
- } // Level-Loop End
-
-
-
+ } // GridLevel-Loop End
//////////////////////////////////////////
//--- Print Storage
+ //////////////////////////////////////////
int tableWidth = 12;
std::cout << center("Levels ",tableWidth) << " | "
<< center("q1",tableWidth) << " | "
@@ -701,10 +354,7 @@ int main(int argc, char *argv[])
std::cout << std::string(tableWidth*9 + 3*9, '-') << "\n";
log << std::string(tableWidth*9 + 3*9, '-') << "\n";
- log.close();
+ log.close(); //close log-file
std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
-
-
-
}
diff --git a/src/deprecated_code/30-8-22_TestingFiles-Comments/Cell-Problem-New.cc b/src/deprecated_code/30-8-22_TestingFiles-Comments/Cell-Problem-New.cc
new file mode 100644
index 0000000000000000000000000000000000000000..22fe2eaf1c788753a83a9618d7f7cf7cc368adbf
--- /dev/null
+++ b/src/deprecated_code/30-8-22_TestingFiles-Comments/Cell-Problem-New.cc
@@ -0,0 +1,798 @@
+#include <config.h>
+#include <array>
+#include <vector>
+#include <fstream>
+
+#include <iostream>
+#include <dune/common/indices.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/math.hh>
+
+
+#include <dune/geometry/quadraturerules.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+// #include <dune/grid/utility/structuredgridfactory.hh> //TEST
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/multitypeblockmatrix.hh>
+#include <dune/istl/multitypeblockvector.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/spqr.hh>
+#include <dune/istl/preconditioners.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/backends/istlvectorbackend.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+#include <dune/functions/functionspacebases/compositebasis.hh>
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/periodicbasis.hh>
+#include <dune/functions/functionspacebases/subspacebasis.hh>
+#include <dune/functions/functionspacebases/boundarydofs.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/microstructure/prestrain_material_geometry.hh>
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/vtk_filler.hh> //TEST
+#include <dune/microstructure/CorrectorComputer.hh>
+#include <dune/microstructure/EffectiveQuantitiesComputer.hh>
+#include <dune/microstructure/prestrainedMaterial.hh>
+
+#include <dune/solvers/solvers/umfpacksolver.hh> //TEST
+#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
+
+// #include <dune/fufem/discretizationerror.hh>
+#include <dune/fufem/dunepython.hh>
+#include <python2.7/Python.h>
+
+#include <dune/fufem/functions/virtualgridfunction.hh> //TEST
+
+// #include <boost/multiprecision/cpp_dec_float.hpp>
+#include <any>
+#include <variant>
+#include <string>
+#include <iomanip> // needed when working with relative paths e.g. from python-scripts
+
+using namespace Dune;
+using namespace MatrixOperations;
+
+//////////////////////////////////////////////////////////////////////
+// Helper functions for Table-Output
+//////////////////////////////////////////////////////////////////////
+/*! Center-aligns string within a field of width w. Pads with blank spaces
+ to enforce alignment. */
+std::string center(const std::string s, const int w) {
+ std::stringstream ss, spaces;
+ int padding = w - s.size(); // count excess room to pad
+ for(int i=0; i<padding/2; ++i)
+ spaces << " ";
+ ss << spaces.str() << s << spaces.str(); // format with padding
+ if(padding>0 && padding%2!=0) // if odd #, add 1 space
+ ss << " ";
+ return ss.str();
+}
+
+/* Convert double to string with specified number of places after the decimal
+ and left padding. */
+template<class type>
+std::string prd(const type x, const int decDigits, const int width) {
+ std::stringstream ss;
+// ss << std::fixed << std::right;
+ ss << std::scientific << std::right; // Use scientific Output!
+ ss.fill(' '); // fill space around displayed #
+ ss.width(width); // set width around displayed #
+ ss.precision(decDigits); // set # places after decimal
+ ss << x;
+ return ss.str();
+}
+
+//////////////////////////////////////////////////
+// Infrastructure for handling periodicity
+//////////////////////////////////////////////////
+// Check whether two points are equal on R/Z x R/Z x R
+auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
+ {
+ return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
+ and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
+ and (FloatCmp::eq(x[2],y[2]))
+ );
+ };
+
+
+
+// // a function:
+// int half(int x, int y) {return x/2+y/2;}
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[])
+{
+ MPIHelper::instance(argc, argv);
+
+ Dune::Timer globalTimer;
+
+ 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");
+
+ //--- setup Log-File
+ std::fstream log;
+ log.open(outputPath + "/output.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;
+
+
+ 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);
+
+
+
+
+ ///////////////////////////////////
+ // Get Parameters/Data
+ ///////////////////////////////////
+ double gamma = parameterSet.get<double>("gamma",1.0); // ratio dimension reduction to homogenization
+ double alpha = parameterSet.get<double>("alpha", 2.0);
+ double theta = parameterSet.get<double>("theta",1.0/4.0);
+ ///////////////////////////////////
+ // Get Material Parameters
+ ///////////////////////////////////
+ std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
+ log << "material_prestrain used: "<< imp << std::endl;
+ double beta = parameterSet.get<double>("beta",2.0);
+ double mu1 = parameterSet.get<double>("mu1",1.0);;
+ double mu2 = beta*mu1;
+ double lambda1 = parameterSet.get<double>("lambda1",0.0);;
+ double lambda2 = beta*lambda1;
+
+
+ if(imp == "material_neukamm")
+ {
+ std::cout <<"mu: " << parameterSet.get<std::array<double,3>>("mu", {1.0,3.0,2.0}) << std::endl;
+ std::cout <<"lambda: " << parameterSet.get<std::array<double,3>>("lambda", {1.0,3.0,2.0}) << std::endl;
+ }
+ else
+ {
+ std::cout <<"mu: " << parameterSet.get<double>("mu1",1.0) << std::endl;
+ std::cout <<"lambda: " << parameterSet.get<double>("lambda1",0.0) << std::endl;
+ }
+
+ ///////////////////////////////////
+ // Get Prestrain/Parameters
+ ///////////////////////////////////
+ auto prestrainImp = PrestrainImp<dim>();
+ auto B_Term = prestrainImp.getPrestrain(parameterSet);
+
+
+
+ log << "----- Input Parameters -----: " << std::endl;
+ log << "alpha: " << alpha << std::endl;
+ log << "gamma: " << gamma << std::endl;
+ log << "theta: " << theta << std::endl;
+ log << "beta: " << beta << std::endl;
+ log << "material parameters: " << std::endl;
+ log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
+ log << "lambda1: " << lambda1 <<"\nlambda2: " << lambda2 << std::endl;
+ log << "----------------------------: " << std::endl;
+
+ ///////////////////////////////////
+ // Generate the grid
+ ///////////////////////////////////
+ // --- Corrector Problem Domain (-1/2,1/2)^3:
+ 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});
+ int levelCounter = 0;
+
+
+ ///////////////////////////////////
+ // 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_Quantities;
+
+ // for(const size_t &level : numLevels) // explixite Angabe der levels.. {2,4}
+ for(size_t level = numLevels[0] ; level <= numLevels[1]; level++) // levels von bis.. [2,4]
+ {
+ std::cout << " ----------------------------------" << std::endl;
+ std::cout << "GridLevel: " << level << std::endl;
+ std::cout << " ----------------------------------" << std::endl;
+
+ 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;
+ log << "Number of Grid-Elements in each direction: " << nElements << std::endl;
+
+ using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+ CellGridType grid_CE(lower,upper,nElements);
+ using GridView = CellGridType::LeafGridView;
+ const GridView gridView_CE = grid_CE.leafGridView();
+ if(print_debug)
+ std::cout << "Host grid has " << gridView_CE.size(dim) << " vertices." << std::endl;
+
+ // //not needed
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ // using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ // using VectorCT = BlockVector<FieldVector<double,1> >;
+ // using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
+
+ ///////////////////////////////////
+ // Create Lambda-Functions for material Parameters depending on microstructure
+ ///////////////////////////////////
+ auto materialImp = IsotropicMaterialImp<dim>();
+ auto muTerm = materialImp.getMu(parameterSet);
+ auto lambdaTerm = materialImp.getLambda(parameterSet);
+
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
+ auto muLocal = localFunction(muGridF);
+ auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
+ auto lambdaLocal = localFunction(lambdaGridF);
+
+
+ //--- Choose a finite element space for Cell Problem
+ 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.)
+ for (const auto& v1 : vertices(gridView_CE))
+ for (const auto& v2 : vertices(gridView_CE))
+ if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
+ {
+ periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
+ }
+
+ //--- setup first order periodic Lagrange-Basis
+ auto Basis_CE = makeBasis(
+ gridView_CE,
+ power<dim>( // eig dimworld?!?!
+ Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
+ flatLexicographic()
+ //blockedInterleaved() // Not Implemented
+ ));
+ if(print_debug)
+ std::cout << "power<periodic> basis has " << Basis_CE.dimension() << " degrees of freedom" << std::endl;
+
+
+
+
+
+ //TEST
+ //Read from Parset...
+ // int Phases = parameterSet.get<int>("Phases", 3);
+
+
+ // std::string materialFunctionName = parameterSet.get<std::string>("materialFunction", "material");
+ // Python::Module module = Python::import(materialFunctionName);
+
+
+
+ // auto indicatorFunction = Python::make_function<double>(module.get("f"));
+ // Func2Tensor indicatorFunction = Python::make_function<double>(module.get("f"));
+ // auto materialFunction_ = Python::make_function<double>(module.get("f"));
+ // auto materialFunction_ = Python::make_function<double>(module.get("f"));
+ // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
+
+ // int Phases;
+ // module.get("Phases").toC<int>(Phases);
+ // std::cout << "Number of Phases used:" << Phases << std::endl;
+
+
+ // std::cout << typeid(mu_).name() << '\n';
+
+
+ //---- Get mu/lambda values (for isotropic material) from Material-file
+ // FieldVector<double,3> mu_(0);
+ // module.get("mu_").toC<FieldVector<double,3>>(mu_);
+ // printvector(std::cout, mu_ , "mu_", "--");
+ // FieldVector<double,3> lambda_(0);
+ // module.get("lambda_").toC<FieldVector<double,3>>(lambda_);
+ // printvector(std::cout, lambda_ , "lambda_", "--");
+
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // TESTING PRESTRAINEDMATERIAL.HH:
+ using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
+
+ auto material_ = prestrainedMaterial(gridView_CE,parameterSet);
+ // Func2Tensor elasticityTensor = material_.getElasticityTensor();
+ // auto elasticityTensor = material_.getElasticityTensor();
+ // Func2TensorParam elasticityTensor_ = *material_.getElasticityTensor();
+ auto elasticityTensor_ = material_.getElasticityTensor();
+
+ // Func2TensorParam TestTensor = Python::make_function<MatrixRT>(module.get("H"));
+
+ // std::cout << "decltype(elasticityTensor_) " << decltype(elasticityTensor_) << std::endl;
+
+
+ // std::cout << "typeid(elasticityTensor).name() :" << typeid(elasticityTensor_).name() << '\n';
+ // std::cout << "typeid(TestTensor).name() :" << typeid(TestTensor).name() << '\n';
+
+ using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
+ // std::cout << "Import NOW:" << std::endl;
+ // MatrixFunc symTest = Python::make_function<MatrixRT>(module.get("sym"));
+ // auto indicatorFunction = Python::make_function<double>(module.get("indicatorFunction"));
+
+ //localize..
+ // auto indicatorFunctionGVF = Dune::Functions::makeGridViewFunction(indicatorFunction , Basis_CE.gridView());
+ // GridViewFunction<double(const Domain&), GridView> indicatorFunctionGVF = Dune::Functions::makeGridViewFunction(indicatorFunction , Basis_CE.gridView());
+ // auto localindicatorFunction = localFunction(indicatorFunctionGVF);
+
+
+ // auto indicatorFunctionGVF = material_.getIndicatorFunction();
+ // auto localindicatorFunction = localFunction(indicatorFunctionGVF);
+
+ auto localindicatorFunction = material_.getIndicatorFunction();
+
+ // GridView::Element
+ // auto localindicatorFunction = localFunction(indicatorFunction);
+
+ // std::cout << "typeid(localindicatorFunction).name() :" << typeid(localindicatorFunction).name() << '\n';
+
+ // using MatrixDomainFunc = std::function< MatrixRT(const MatrixRT&,const Domain&)>;
+ // // MatrixFunc elasticityTensor = Python::make_function<MatrixRT>(module.get("L"));
+
+ // auto elasticityTensorGVF = Dune::Functions::makeGridViewFunction(elasticityTensor , Basis_CE.gridView());
+ // auto localElasticityTensor = localFunction(elasticityTensorGVF);
+
+ // Func2Tensor forceTerm = [] (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?
+ // };
+
+ // auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, Basis_CE.gridView());
+ // auto loadFunctional = localFunction(loadGVF);
+
+ MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+ // auto xu = symTest(G1_);
+ // std::cout << "TEST NOW:" << std::endl;
+ // printmatrix(std::cout, symTest(G1_), "symTest(G1_)", "--");
+
+ // auto TestTensorGVF = Dune::Functions::makeGridViewFunction(TestTensor , Basis_CE.gridView());
+ // auto localTestTensor = localFunction(TestTensorGVF );
+
+
+ // printmatrix(std::cout, elasticityTensor(G1_), "elasticityTensor(G1_)", "--");
+ // auto temp = elasticityTensor(G1_);
+
+
+
+ // for (const auto& element : elements(Basis_CE.gridView()))
+ // {
+ // localindicatorFunction.bind(element);
+ // int orderQR = 2;
+ // const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ // for (const auto& quadPoint : quad)
+ // {
+ // const auto& quadPos = quadPoint.position();
+
+ // // std::cout << "localindicatorFunction(quadPos): " << localindicatorFunction(quadPos) << std::endl;
+
+
+
+ // // std::cout << "quadPos : " << quadPos << std::endl;
+ // // auto temp = TestTensor(G1_, element.geometry().global(quadPos));
+ // // auto temp2 = elasticityTensor_(G1_, element.geometry().global(quadPos));
+ // // std::cout << "material_.applyElasticityTensor:" << std::endl;
+ // // auto tmp3 = material_.applyElasticityTensor(G1_, element.geometry().global(quadPos));
+ // // auto tmp3 = material_.applyElasticityTensor(G1_, quadPos);
+ // // printmatrix(std::cout, tmp3, "tmp3", "--");
+ // }
+ // }
+
+
+
+ // for (auto&& vertex : vertices(gridView_CE))
+ // {
+ // std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
+ // auto tmp = vertex.geometry().corner(0);
+ // auto temp = elasticityTensor(tmp);
+ // // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
+ // // printmatrix(std::cout, localElasticityTensor(G1_,tmp), "localElasticityTensor(vertex.geometry().corner(0))", "--");
+ // }
+
+
+
+ // std::function<int(int,int)> fn1 = half;
+ // std::cout << "fn1(60,20): " << fn1(60,20) << '\n';
+
+
+ // std::cout << typeid(elasticityTensorGVF).name() << '\n';
+ // std::cout << typeid(localElasticityTensor).name() << '\n';
+
+
+ // ParameterTree parameterSet_2;
+ // ParameterTreeParser::readINITree(geometryFunctionPath + "/"+ materialFunctionName + ".py", parameterSet_2);
+
+ // auto lu = parameterSet_2.get<FieldVector<double,3>>("lu", {1.0,3.0,2.0});
+ // std::cout <<"lu[1]: " << lu[1]<< std::endl;
+ // std::cout <<"lu: " << parameterSet_2.get<std::array<double,3>>("lu", {1.0,3.0,2.0}) << std::endl;
+
+ // auto mU_ = module.evaluate(parameterSet_2.get<std::string>("lu", "[1,2,3]"));
+ // std::cout << "typeid(mU_).name()" << typeid(mU_.operator()()).name() << '\n';
+
+
+ // for (auto&& vertex : vertices(gridView_CE))
+ // {
+ // std::cout << "vertex.geometry().corner(0):" << vertex.geometry().corner(0)<< std::endl;
+ // // std::cout << "materialFunction_(vertex.geometry().corner(0))", materialFunction_(vertex.geometry().corner(0)) << std::endl;
+ // printvector(std::cout, materialFunction_(vertex.geometry().corner(0)), "materialFunction_(vertex.geometry().corner(0))", "--");
+ // }
+ // std::cout << "materialFunction_({0.0,0.0,0.0})", materialFunction_({0.0,0.0,0.0}) << std::endl;
+
+
+
+ // --------------------------------------------------------------
+
+ //TODO// Phasen anhand von Mu bestimmen?
+ //TODO: DUNE_THROW(Exception, "Inconsistent choice of interpolation method"); if number of Phases != mu/lambda parameters
+
+ //FÜR L GARNICHT NÖTIG DENN RÜCKGABETYPE IS IMMER MATRIXRT!?!:
+ // BEi materialfunction (isotopic) reicht auch FieldVector<double,2> für lambda/mu
+
+ // switch (Phases)
+ // {
+ // case 1: //homogeneous material
+ // {
+ // std::cout << "Phase - 1" << std::endl;
+ // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
+ // break;
+ // }
+ // case 2:
+ // {
+ // std::cout << "Phase - 1" << std::endl;
+ // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
+ // break;
+ // }
+ // case 3:
+ // {
+ // std::cout << "Phase - 3" << std::endl;
+ // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
+ // break;
+ // }
+ // }
+
+
+ // switch (Phases)
+ // {
+ // case 1: //homogeneous material
+ // {
+ // std::cout << "Phases - 1" << std::endl;
+ // std::array<double,1> mu_ = parameterSet.get<std::array<double,1>>("mu", {1.0});
+ // Python::Module module = Python::import(materialFunction);
+ // auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
+ // auto muTerm = [mu_] (const Domain& x) {return mu_;};
+ // break;
+ // }
+ // case 2:
+ // {
+ // std::cout << "Phases - 2" << std::endl;
+ // std::array<double,2> mu_ = parameterSet.get<std::array<double,2>>("mu", {1.0,3.0});
+ // Python::Module module = Python::import(materialFunction);
+ // auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
+ // auto muTerm = [mu_,indicatorFunction] (const Domain& x)
+ // {
+ // if (indicatorFunction(x) == 1)
+ // return mu_[0];
+ // else
+ // return mu_[1];
+ // };
+ // break;
+ // }
+ // case 3:
+ // {
+ // std::cout << "Phases - 3" << std::endl;
+ // std::array<double,3> mu_ = parameterSet.get<std::array<double,3>>("mu", {1.0,3.0, 5.0});
+ // Python::Module module = Python::import(materialFunction);
+ // auto indicatorFunction = Python::make_function<double>(module.get("f")); // get indicator function
+ // auto muTerm = [mu_,indicatorFunction] (const Domain& x)
+ // {
+ // if (indicatorFunction(x) == 1)
+ // return mu_[0];
+ // else if (indicatorFunction(x) == 2)
+ // return mu_[1];
+ // else
+ // return mu_[2];
+ // };
+ // break;
+ // }
+ // }
+
+
+
+
+ //TEST
+// std::cout << "Test crossSectionDirectionScaling:" << std::endl;
+/*
+ MatrixRT T {{1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}};
+ printmatrix(std::cout, T, "Matrix T", "--");
+
+ auto ST = crossSectionDirectionScaling((1.0/5.0),T);
+ printmatrix(std::cout, ST, "scaled Matrix T", "--");*/
+
+ //TEST
+// auto QuadraticForm = [] (const double mu, const double lambda, const MatrixRT& M) {
+//
+// return lambda*std::pow(trace(M),2) + 2*mu*pow(norm(sym(M)),2);
+// };
+
+
+
+
+ //------------------------------------------------------------------------------------------------
+ //--- compute Correctors
+ // auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
+ auto correctorComputer = CorrectorComputer(Basis_CE, material_, muTerm, lambdaTerm, gamma, log, parameterSet);
+ correctorComputer.solve();
+
+
+
+//////////////////////////////////////////////////
+
+
+ //--- check Correctors (options):
+ if(parameterSet.get<bool>("write_L2Error", false))
+ correctorComputer.computeNorms();
+ if(parameterSet.get<bool>("write_VTK", false))
+ correctorComputer.writeCorrectorsVTK(level);
+ //--- additional Test: check orthogonality (75) from paper:
+ if(parameterSet.get<bool>("write_checkOrthogonality", false))
+ correctorComputer.check_Orthogonality();
+ //--- Check symmetry of stiffness matrix
+ if(print_debug)
+ correctorComputer.checkSymmetry();
+
+ // auto B_Term = material_.getPrestrainFunction();
+
+
+ //--- compute effective quantities
+ auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term,material_);
+ effectiveQuantitiesComputer.computeEffectiveQuantities();
+
+
+
+
+ //--- Test:: Compute Qeff without using the orthogonality (75)...
+ // only really makes a difference whenever the orthogonality is not satisfied!
+ // std::cout << "----------computeFullQ-----------"<< std::endl; //TEST
+ // effectiveQuantitiesComputer.computeFullQ();
+
+ //--- get effective quantities
+ auto Qeff = effectiveQuantitiesComputer.getQeff();
+ auto Beff = effectiveQuantitiesComputer.getBeff();
+ printmatrix(std::cout, Qeff, "Matrix Qeff", "--");
+ printvector(std::cout, Beff, "Beff", "--");
+
+ //--- write effective quantities to matlab folder (for symbolic minimization)
+ if(parameterSet.get<bool>("write_toMATLAB", false))
+ effectiveQuantitiesComputer.writeToMatlab(outputPath);
+
+ std::cout.precision(10);
+ std::cout<< "q1 : " << Qeff[0][0] << std::endl;
+ std::cout<< "q2 : " << Qeff[1][1] << std::endl;
+ std::cout<< "q3 : " << std::fixed << Qeff[2][2] << std::endl;
+ std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Qeff[0][1] << std::endl;
+ // -------------------------------------------
+
+
+ //TEST
+ // 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?
+ // };
+
+ // double energy = effectiveQuantitiesComputer.energySP(x3G_1,x3G_1);
+ // std::cout << "energy:" << energy << std::endl;
+
+ Storage_Quantities.push_back(Qeff[0][0] );
+ 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[1][2] );
+ Storage_Quantities.push_back(Beff[0]);
+ Storage_Quantities.push_back(Beff[1]);
+ Storage_Quantities.push_back(Beff[2]);
+
+ log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
+ log << "q1=" << Qeff[0][0] << std::endl;
+ log << "q2=" << Qeff[1][1] << std::endl;
+ log << "q3=" << Qeff[2][2] << std::endl;
+ log << "q12=" << Qeff[0][1] << 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;
+ log << "b2=" << Beff[1] << std::endl;
+ log << "b3=" << Beff[2] << std::endl;
+ log << "mu_gamma=" << Qeff[2][2] << std::endl; // added for Python-Script
+
+
+
+
+ if (write_materialFunctions)
+ {
+ using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+// VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
+// VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{40,40,40});
+ VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
+
+ using GridViewVTK = VTKGridType::LeafGridView;
+ const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+
+ auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+ auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+ std::vector<double> mu_CoeffP0;
+ Functions::interpolate(scalarP0FeBasis, mu_CoeffP0, muTerm);
+ auto mu_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, mu_CoeffP0);
+
+ std::vector<double> mu_CoeffP1;
+ Functions::interpolate(scalarP1FeBasis, mu_CoeffP1, muTerm);
+ auto mu_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, mu_CoeffP1);
+
+
+ std::vector<double> lambda_CoeffP0;
+ Functions::interpolate(scalarP0FeBasis, lambda_CoeffP0, lambdaTerm);
+ auto lambda_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, lambda_CoeffP0);
+
+ std::vector<double> lambda_CoeffP1;
+ Functions::interpolate(scalarP1FeBasis, lambda_CoeffP1, lambdaTerm);
+ auto lambda_DGBF_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, lambda_CoeffP1);
+
+ VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
+
+ MaterialVtkWriter.addVertexData(
+ mu_DGBF_P1,
+ VTK::FieldInfo("mu_P1", VTK::FieldInfo::Type::scalar, 1));
+ MaterialVtkWriter.addCellData(
+ mu_DGBF_P0,
+ VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));
+ MaterialVtkWriter.addVertexData(
+ lambda_DGBF_P1,
+ VTK::FieldInfo("lambda_P1", VTK::FieldInfo::Type::scalar, 1));
+ MaterialVtkWriter.addCellData(
+ lambda_DGBF_P0,
+ VTK::FieldInfo("lambda_P0", 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;
+
+ }
+
+
+// if (write_prestrainFunctions)
+// {
+// using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+// // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
+// // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{40,40,40});
+// VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
+// using GridViewVTK = VTKGridType::LeafGridView;
+// const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
+
+// FTKfillerContainer<dim> VTKFiller;
+// VTKFiller.vtkPrestrainNorm(gridView_VTK, B_Term, "PrestrainBNorm");
+
+// // WORKS Too
+// VTKFiller.vtkProblemCell(gridView_VTK, B_Term, muLocal,"VTKProblemCell");;
+
+
+// // TEST
+// auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
+// auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
+
+// std::vector<double> B_CoeffP0;
+// Functions::interpolate(scalarP0FeBasis, B_CoeffP0, B_Term);
+// auto B_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, B_CoeffP0);
+
+
+
+// VTKWriter<GridView> PrestrainVtkWriter(gridView_VTK);
+
+// PrestrainVtkWriter.addCellData(
+// B_DGBF_P0,
+// VTK::FieldInfo("B_P0", VTK::FieldInfo::Type::scalar, 1));
+
+// PrestrainVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
+// std::cout << "wrote data to file:" + outputPath +"/PrestrainFunctions-level" + std::to_string(level) << std::endl;
+
+// }
+
+
+
+
+ levelCounter++;
+ } // Level-Loop End
+
+
+
+
+ //////////////////////////////////////////
+ //--- Print Storage
+ int tableWidth = 12;
+ std::cout << center("Levels ",tableWidth) << " | "
+ << center("q1",tableWidth) << " | "
+ << center("q2",tableWidth) << " | "
+ << center("q3",tableWidth) << " | "
+ << center("q12",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";
+ log << center("Levels ",tableWidth) << " | "
+ << center("q1",tableWidth) << " | "
+ << center("q2",tableWidth) << " | "
+ << center("q3",tableWidth) << " | "
+ << center("q12",tableWidth) << " | "
+ << center("q23",tableWidth) << " | "
+ << center("b1",tableWidth) << " | "
+ << center("b2",tableWidth) << " | "
+ << center("b3",tableWidth) << " | " << "\n";
+ log << std::string(tableWidth*9 + 3*9, '-') << "\n";
+
+ int StorageCount2 = 0;
+ for(auto& v: Storage_Quantities)
+ {
+ 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 )
+ {
+ std::cout << std::endl;
+ log << std::endl;
+ }
+ }
+ std::cout << std::string(tableWidth*9 + 3*9, '-') << "\n";
+ log << std::string(tableWidth*9 + 3*9, '-') << "\n";
+
+ log.close();
+
+ std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
+
+
+
+}
diff --git a/src/Cell-Problem.cc b/src/deprecated_code/30-8-22_TestingFiles-Comments/Cell-Problem.cc
similarity index 100%
rename from src/Cell-Problem.cc
rename to src/deprecated_code/30-8-22_TestingFiles-Comments/Cell-Problem.cc
diff --git a/src/deprecated_code/30-8-22_TestingFiles-Comments/CorrectorComputer.hh b/src/deprecated_code/30-8-22_TestingFiles-Comments/CorrectorComputer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..762125683373d21057d427b7c73397f7168da99e
--- /dev/null
+++ b/src/deprecated_code/30-8-22_TestingFiles-Comments/CorrectorComputer.hh
@@ -0,0 +1,1445 @@
+#ifndef DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
+#define DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
+
+#include <dune/common/parametertree.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/microstructure/matrix_operations.hh>
+
+#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
+#include <dune/solvers/solvers/umfpacksolver.hh>
+
+using namespace Dune;
+// using namespace Functions;
+using namespace MatrixOperations;
+
+using std::shared_ptr;
+using std::make_shared;
+using std::fstream;
+
+
+template <class Basis, class Material> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
+class CorrectorComputer {
+
+public:
+ static const int dimworld = 3; //GridView::dimensionworld;
+ static const int dim = Basis::GridView::dimension; //const int dim = Domain::dimension;
+
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+
+ using ScalarRT = FieldVector< double, 1>;
+ using VectorRT = FieldVector< double, dimworld>;
+ using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+
+ 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> >;
+ using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
+
+ using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
+
+protected:
+//private:
+ const Basis& basis_;
+
+
+ const Material& material_;
+
+ fstream& log_; // Output-log
+ const ParameterTree& parameterSet_;
+
+ const FuncScalar mu_;
+ const FuncScalar lambda_;
+ double gamma_;
+
+ MatrixCT stiffnessMatrix_;
+ VectorCT load_alpha1_,load_alpha2_,load_alpha3_; //right-hand side(load) vectors
+
+ 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_};
+
+ // --- Offset between basis indices
+ const int phiOffset_;
+
+public:
+ ///////////////////////////////
+ // constructor
+ ///////////////////////////////
+ CorrectorComputer( const Basis& basis,
+ const Material& material,
+ const FuncScalar& mu,
+ const FuncScalar& lambda,
+ double gamma,
+ std::fstream& log,
+ const ParameterTree& parameterSet)
+ : basis_(basis),
+ material_(material),
+ mu_(mu),
+ lambda_(lambda),
+ gamma_(gamma),
+ log_(log),
+ parameterSet_(parameterSet),
+ phiOffset_(basis.size())
+ {
+
+ assemble();
+
+ // if (parameterSet.get<bool>("stiffnessmatrix_cellproblem_to_csv"))
+ // csvSystemMatrix();
+ // if (parameterSet.get<bool>("rhs_cellproblem_to_csv"))
+ // csvRHSs();
+ // if (parameterSet.get<bool>("rhs_cellproblem_to_vtk"))
+ // vtkLoads();
+ }
+
+
+ // -----------------------------------------------------------------
+ // --- Assemble Corrector problems
+ void assemble()
+ {
+ Dune::Timer StiffnessTimer;
+ assembleCellStiffness(stiffnessMatrix_);
+ std::cout << "Stiffness assembly Timer: " << StiffnessTimer.elapsed() << std::endl;
+
+ assembleCellLoad(load_alpha1_ ,x3G_1_);
+ assembleCellLoad(load_alpha2_ ,x3G_2_);
+ assembleCellLoad(load_alpha3_ ,x3G_3_);
+ };
+
+
+
+ ///////////////////////////////
+ // getter
+ ///////////////////////////////
+ const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}
+
+ ParameterTree getParameterSet() const {return parameterSet_;}
+
+ fstream* getLog(){return &log_;}
+
+ double getGamma(){return gamma_;}
+
+ 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_);}
+
+ shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
+
+ // --- 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_);}
+
+
+
+
+ // Get the occupation pattern of the stiffness matrix
+ void getOccupationPattern(MatrixIndexSet& nb)
+ {
+ // OccupationPattern:
+ // | phi*phi m*phi |
+ // | phi *m m*m |
+ auto localView = basis_.localView();
+ const int phiOffset = basis_.dimension();
+
+ nb.resize(basis_.size()+3, basis_.size()+3);
+
+ for (const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ for (size_t i=0; i<localView.size(); i++)
+ {
+ // The global index of the i-th vertex of the element
+ auto row = localView.index(i);
+ for (size_t j=0; j<localView.size(); j++ )
+ {
+ // The global index of the j-th vertex of the element
+ auto col = localView.index(j);
+ nb.add(row[0],col[0]); // nun würde auch nb.add(row,col) gehen..
+ }
+ }
+ for (size_t i=0; i<localView.size(); i++)
+ {
+ auto row = localView.index(i);
+
+ for (size_t j=0; j<3; j++ )
+ {
+ nb.add(row,phiOffset+j); // m*phi part of StiffnessMatrix
+ nb.add(phiOffset+j,row); // phi*m part of StiffnessMatrix
+ }
+ }
+ for (size_t i=0; i<3; i++ )
+ for (size_t j=0; j<3; j++ )
+ {
+ nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
+ }
+ }
+ //////////////////////////////////////////////////////////////////
+ // setOneBaseFunctionToZero
+ //////////////////////////////////////////////////////////////////
+ if(parameterSet_.get<bool>("set_oneBasisFunction_Zero ", true)){
+ FieldVector<int,3> row;
+ unsigned int arbitraryLeafIndex = parameterSet_.get<unsigned int>("arbitraryLeafIndex", 0);
+ unsigned int arbitraryElementNumber = parameterSet_.get<unsigned int>("arbitraryElementNumber", 0);
+
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
+ const auto nSf = localFiniteElement.localBasis().size();
+ assert(arbitraryLeafIndex < nSf );
+ assert(arbitraryElementNumber < basis_.gridView().size(0)); // "arbitraryElementNumber larger than total Number of Elements"
+
+ //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+ row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
+
+ for (int k = 0; k<3; k++)
+ nb.add(row[k],row[k]);
+ }
+ std::cout << "finished occupation pattern\n";
+ }
+
+
+ template<class localFunction1, class localFunction2>
+ void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
+ ElementMatrixCT& elementMatrix,
+ const localFunction1& mu,
+ const localFunction2& lambda
+ )
+ {
+ // Local StiffnessMatrix of the form:
+ // | phi*phi m*phi |
+ // | phi *m m*m |
+ auto element = localView.element();
+ auto geometry = element.geometry();
+ // using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+
+ elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
+ elementMatrix = 0;
+
+
+ auto elasticityTensor = material_.getElasticityTensor();
+ // 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
+ const int localPhiOffset = localView.size();
+
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+ // std::cout << "localView.size(): " << localView.size() << std::endl;
+ // std::cout << "nSf : " << nSf << std::endl;
+
+ ///////////////////////////////////////////////
+ // Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
+ //////////////////////////////////////////////
+ MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{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 > basisContainer = {G_1, G_2, G_3};
+
+ // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ // int orderQR = 0;
+ // int orderQR = 1;
+ // int orderQR = 2;
+ // int orderQR = 3;
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+ // double elementContribution = 0.0;
+
+ // std::cout << "Print QuadratureOrder:" << orderQR << std::endl; //TEST`
+
+ int QPcounter= 0;
+ for (const auto& quadPoint : quad)
+ {
+ QPcounter++;
+ const auto& quadPos = quadPoint.position();
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+ const auto integrationElement = geometry.integrationElement(quadPos);
+
+ std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,
+ referenceGradients);
+ // Compute the shape function gradients on the grid element
+ std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+
+ for (size_t i=0; i<gradients.size(); i++)
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+
+ for (size_t l=0; l< dimworld; l++)
+ for (size_t j=0; j < nSf; j++ )
+ {
+ size_t row = localView.tree().child(l).localIndex(j);
+ // (scaled) Deformation gradient of the test basis function
+ MatrixRT defGradientV(0);
+ defGradientV[l][0] = gradients[j][0]; // Y
+ defGradientV[l][1] = gradients[j][1]; //X2
+ // defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
+ defGradientV[l][2] = gradients[j][2]; //X3
+
+ defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
+ // "phi*phi"-part
+ for (size_t k=0; k < dimworld; k++)
+ for (size_t i=0; i < nSf; i++)
+ {
+ // (scaled) Deformation gradient of the ansatz basis function
+ MatrixRT defGradientU(0);
+ defGradientU[k][0] = gradients[i][0]; // Y
+ defGradientU[k][1] = gradients[i][1]; //X2
+ // defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
+ defGradientU[k][2] = gradients[i][2]; //X3
+ // printmatrix(std::cout, defGradientU , "defGradientU", "--");
+ defGradientU = crossSectionDirectionScaling((1.0/gamma_),defGradientU);
+
+
+ // auto etmp = material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos));
+ // auto etmp = elasticityTensor(defGradientU,element.geometry().global(quadPos));
+ // auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos);
+ // printmatrix(std::cout, etmp, "etmp", "--");
+ // double energyDensity= scalarProduct(etmp,defGradientV);
+ // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),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_.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));
+ // 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);
+
+
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
+ // double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
+
+ size_t col = localView.tree().child(k).localIndex(i);
+
+ elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
+ }
+
+ // "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_.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;
+ // 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(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);
+ // double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST
+ auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
+ elementMatrix[row][localPhiOffset+m] += value;
+ elementMatrix[localPhiOffset+m][row] += value;
+ }
+
+ }
+ // "m*m"-part
+ for(size_t m=0; m<3; m++) //TODO ist symmetric.. reicht die hälfte zu berechnen!!!
+ for(size_t n=0; n<3; n++)
+ {
+
+ // std::cout << "QPcounter: " << QPcounter << std::endl;
+ // std::cout << "m= " << m << " n = " << n << std::endl;
+ // printmatrix(std::cout, basisContainer[m] , "basisContainer[m]", "--");
+ // printmatrix(std::cout, basisContainer[n] , "basisContainer[n]", "--");
+ // std::cout << "integrationElement:" << integrationElement << std::endl;
+ // std::cout << "quadPoint.weight(): "<< quadPoint.weight() << std::endl;
+ // std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
+ // std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
+
+ // 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_.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]));
+ // 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]);
+ // double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST
+ elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
+
+ // std::cout << "energyDensityGG:" << energyDensityGG<< std::endl;
+ // std::cout << "product " << energyDensityGG * quadPoint.weight() * integrationElement << std::endl;
+ // printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
+ }
+ }
+ // std::cout << "Number of QuadPoints:" << QPcounter << std::endl;
+ // printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
+ }
+
+
+
+ // Compute the source term for a single element
+ // < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
+ template<class LocalFunction1, class LocalFunction2, class Vector, class LocalForce>
+ void computeElementLoadVector( const typename Basis::LocalView& localView,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
+ Vector& elementRhs,
+ const LocalForce& forceTerm
+ )
+ {
+ // | f*phi|
+ // | --- |
+ // | f*m |
+ // using Element = typename LocalView::Element;
+ const auto element = localView.element();
+ const auto geometry = element.geometry();
+ // constexpr int dim = Element::dimension;
+ // constexpr int dimworld = dim;
+
+ // using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+
+ auto elasticityTensor = material_.getElasticityTensor();
+ // 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();
+ const auto nSf = localFiniteElement.localBasis().size();
+
+ elementRhs.resize(localView.size() +3);
+ elementRhs = 0;
+
+ // LocalBasis-Offset
+ const int localPhiOffset = localView.size();
+
+ ///////////////////////////////////////////////
+ // Basis for R_sym^{2x2}
+ //////////////////////////////////////////////
+ MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{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 > basisContainer = {G_1, G_2, G_3};
+
+ // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
+ // int orderQR = 0;
+ // int orderQR = 1;
+ // int orderQR = 2;
+ // int orderQR = 3;
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ // std::cout << "Quad-Rule order used: " << orderQR << std::endl;
+
+ for (const auto& quadPoint : quad)
+ {
+ const FieldVector<double,dim>& quadPos = quadPoint.position();
+ const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+ std::vector<FieldMatrix<double,1,dim> > referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+
+ std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());
+ for (size_t i=0; i< gradients.size(); i++)
+ jacobian.mv(referenceGradients[i][0], gradients[i]);
+
+ //TEST
+ // std::cout << "forceTerm(element.geometry().global(quadPos)):" << std::endl;
+ // std::cout << forceTerm(element.geometry().global(quadPos)) << std::endl;
+ // std::cout << "forceTerm(quadPos)" << std::endl;
+ // std::cout << forceTerm(quadPos) << std::endl;
+ //
+ // //TEST QUadrature
+ // std::cout << "quadPos:" << quadPos << std::endl;
+ // std::cout << "element.geometry().global(quadPos):" << element.geometry().global(quadPos) << std::endl;
+ // // //
+ // // //
+ // std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
+ // std::cout << "integrationElement(quadPos):" << integrationElement << std::endl;
+ // std::cout << "mu(quadPos) :" << mu(quadPos) << std::endl;
+ // std::cout << "lambda(quadPos) :" << lambda(quadPos) << std::endl;
+
+
+ // "f*phi"-part
+ for (size_t i=0; i < nSf; i++)
+ for (size_t k=0; k < dimworld; k++)
+ {
+ // Deformation gradient of the ansatz basis function
+ MatrixRT defGradientV(0);
+ defGradientV[k][0] = gradients[i][0]; // Y
+ defGradientV[k][1] = gradients[i][1]; // X2
+ // defGradientV[k][2] = (1.0/gamma_)*gradients[i][2]; //
+ defGradientV[k][2] = gradients[i][2]; // X3
+
+ defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
+
+ // double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(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));
+ // 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);
+
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV );
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
+
+ size_t row = localView.tree().child(k).localIndex(i);
+ elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
+ }
+
+ // "f*m"-part
+ for (size_t m=0; m<3; m++)
+ {
+ // double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(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]));
+ // 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
+ // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
+ elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
+ // std::cout << "energyDensityfG * quadPoint.weight() * integrationElement: " << energyDensityfG * quadPoint.weight() * integrationElement << std::endl;
+ }
+ }
+ }
+
+
+ void assembleCellStiffness(MatrixCT& matrix)
+ {
+ std::cout << "assemble Stiffness-Matrix begins." << std::endl;
+
+ MatrixIndexSet occupationPattern;
+ getOccupationPattern(occupationPattern);
+ occupationPattern.exportIdx(matrix);
+ matrix = 0;
+
+ auto localView = basis_.localView();
+ const int phiOffset = basis_.dimension();
+
+ auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ auto muLocal = localFunction(muGridF);
+ auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ auto lambdaLocal = localFunction(lambdaGridF);
+
+ for (const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ muLocal.bind(element);
+ lambdaLocal.bind(element);
+ const int localPhiOffset = localView.size();
+ // Dune::Timer Time;
+ //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+ Matrix<FieldMatrix<double,1,1> > elementMatrix;
+ computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
+
+ // std::cout << "local assembly time:" << Time.elapsed() << std::endl;
+
+ //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+ //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
+ //std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
+
+ //TEST
+ //Check Element-Stiffness-Symmetry:
+ for (size_t i=0; i<localPhiOffset; i++)
+ for (size_t j=0; j<localPhiOffset; j++ )
+ {
+ if(abs(elementMatrix[i][j] - elementMatrix[j][i]) > 1e-12 )
+ std::cout << "ELEMENT-STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
+ }
+ //////////////////////////////////////////////////////////////////////////////
+ // GLOBAL STIFFNES ASSEMBLY
+ //////////////////////////////////////////////////////////////////////////////
+ for (size_t i=0; i<localPhiOffset; i++)
+ for (size_t j=0; j<localPhiOffset; j++ )
+ {
+ auto row = localView.index(i);
+ auto col = localView.index(j);
+ matrix[row][col] += elementMatrix[i][j];
+ }
+ for (size_t i=0; i<localPhiOffset; i++)
+ for(size_t m=0; m<3; m++)
+ {
+ auto row = localView.index(i);
+ matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m];
+ matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
+ }
+ for (size_t m=0; m<3; m++ )
+ for (size_t n=0; n<3; n++ )
+ matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
+
+ // printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
+ }
+ }
+
+
+ void assembleCellLoad(VectorCT& b,
+ const Func2Tensor& forceTerm
+ )
+ {
+ // std::cout << "assemble load vector." << std::endl;
+ b.resize(basis_.size()+3);
+ b = 0;
+
+ auto localView = basis_.localView();
+ const int phiOffset = basis_.dimension();
+
+ // Transform G_alpha's to GridViewFunctions/LocalFunctions
+ auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis_.gridView());
+ auto loadFunctional = localFunction(loadGVF);
+
+ auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ auto muLocal = localFunction(muGridF);
+ auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ auto lambdaLocal = localFunction(lambdaGridF);
+
+
+ // int counter = 1;
+ for (const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ muLocal.bind(element);
+ lambdaLocal.bind(element);
+ loadFunctional.bind(element);
+
+ const int localPhiOffset = localView.size();
+ // std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+
+ VectorCT elementRhs;
+ // std::cout << "----------------------------------Element : " << counter << std::endl; //TEST
+ // counter++;
+ computeElementLoadVector(localView, muLocal, lambdaLocal, elementRhs, loadFunctional);
+ // computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
+ // printvector(std::cout, elementRhs, "elementRhs", "--");
+ // printvector(std::cout, elementRhs, "elementRhs", "--");
+ //////////////////////////////////////////////////////////////////////////////
+ // GLOBAL LOAD ASSEMBLY
+ //////////////////////////////////////////////////////////////////////////////
+ for (size_t p=0; p<localPhiOffset; p++)
+ {
+ auto row = localView.index(p);
+ b[row] += elementRhs[p];
+ }
+ for (size_t m=0; m<3; m++ )
+ b[phiOffset+m] += elementRhs[localPhiOffset+m];
+ //printvector(std::cout, b, "b", "--");
+ }
+ }
+
+ // -----------------------------------------------------------------
+ // --- Functions for global integral mean equals zero constraint
+ auto arbitraryComponentwiseIndices(const int elementNumber,
+ const int leafIdx
+ )
+ {
+ // (Local Approach -- works for non Lagrangian-Basis too)
+ // Input : elementNumber & localIdx
+ // Output : determine all Component-Indices that correspond to that basis-function
+ auto localView = basis_.localView();
+
+ FieldVector<int,3> row;
+ int elementCounter = 0;
+
+ for (const auto& element : elements(basis_.gridView()))
+ {
+ if(elementCounter == elementNumber) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
+ {
+ localView.bind(element);
+
+ for (int k = 0; k < 3; k++)
+ {
+ auto rowLocal = localView.tree().child(k).localIndex(leafIdx); //Input: LeafIndex! TODO bräuchte hier (Inverse ) Local-to-Leaf Idx Map
+ row[k] = localView.index(rowLocal);
+ // std::cout << "rowLocal:" << rowLocal << std::endl;
+ // std::cout << "row[k]:" << row[k] << std::endl;
+ }
+ }
+ elementCounter++;
+ }
+ return row;
+ }
+
+ void setOneBaseFunctionToZero()
+ {
+ std::cout << "Setting one Basis-function to zero" << std::endl;
+
+ // constexpr int dim = Basis::LocalView::Element::dimension;
+
+ unsigned int arbitraryLeafIndex = parameterSet_.template get<unsigned int>("arbitraryLeafIndex", 0);
+ unsigned int arbitraryElementNumber = parameterSet_.template get<unsigned int>("arbitraryElementNumber", 0);
+ //Determine 3 global indices (one for each component of an arbitrary local FE-function)
+ FieldVector<int,3> row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
+
+ for (int k = 0; k<dim; k++)
+ {
+ load_alpha1_[row[k]] = 0.0;
+ load_alpha2_[row[k]] = 0.0;
+ load_alpha3_[row[k]] = 0.0;
+ auto cIt = stiffnessMatrix_[row[k]].begin();
+ auto cEndIt = stiffnessMatrix_[row[k]].end();
+ for (; cIt!=cEndIt; ++cIt)
+ *cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
+ }
+ }
+
+
+ auto childToIndexMap(const int k)
+ {
+ // Input : child/component
+ // Output : determine all Indices that belong to that component
+ auto localView = basis_.localView();
+
+ std::vector<int> r = { };
+ // for (int n: r)
+ // std::cout << n << ","<< std::endl;
+
+ // Determine global indizes for each component k = 1,2,3.. in order to subtract correct (component of) integral Mean
+ // (global) Indices that correspond to component k = 1,2,3
+ for(const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ const auto& localFiniteElement = localView.tree().child(k).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+
+ for(size_t j=0; j<nSf; j++)
+ {
+ auto Localidx = localView.tree().child(k).localIndex(j); // local indices
+ r.push_back(localView.index(Localidx)); // global indices
+ }
+ }
+ // Delete duplicate elements
+ // first remove consecutive (adjacent) duplicates
+ auto last = std::unique(r.begin(), r.end());
+ r.erase(last, r.end());
+ // sort followed by unique, to remove all duplicates
+ std::sort(r.begin(), r.end());
+ last = std::unique(r.begin(), r.end());
+ r.erase(last, r.end());
+ return r;
+ }
+
+
+ auto integralMean(VectorCT& coeffVector)
+ {
+ auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis_,coeffVector);
+ auto localfun = localFunction(GVFunction);
+
+ auto localView = basis_.localView();
+
+ FieldVector<double,3> r = {0.0, 0.0, 0.0};
+ double area = 0.0;
+
+ // Compute Area integral & integral of FE-function
+ for(const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ localfun.bind(element);
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+ // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; //TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+
+ for(const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = element.geometry().integrationElement(quadPos);
+ area += quadPoint.weight() * integrationElement;
+
+ r += localfun(quadPos) * quadPoint.weight() * integrationElement;
+ }
+ }
+ // std::cout << "Domain-Area: " << area << std::endl;
+ return (1.0/area) * r ;
+ }
+
+
+ auto subtractIntegralMean(VectorCT& coeffVector)
+ {
+ // Substract correct Integral mean from each associated component function
+ auto IM = integralMean(coeffVector);
+
+ for(size_t k=0; k<dim; k++)
+ {
+ //std::cout << "Integral-Mean: " << IM[k] << std::endl;
+ auto idx = childToIndexMap(k);
+ for ( int i : idx)
+ coeffVector[i] -= IM[k];
+ }
+ }
+
+
+ // -----------------------------------------------------------------
+ // --- Solving the Corrector-problem:
+ auto solve()
+ {
+ std::cout << "start corrector solver..." << std::endl;
+ bool set_oneBasisFunction_Zero = parameterSet_.get<bool>("set_oneBasisFunction_Zero", false);
+ bool substract_integralMean = false;
+ if(parameterSet_.get<bool>("set_IntegralZero", false))
+ {
+ set_oneBasisFunction_Zero = true;
+ substract_integralMean = true;
+ }
+ // set one basis-function to zero
+ if(set_oneBasisFunction_Zero)
+ setOneBaseFunctionToZero();
+
+ //TEST: Compute Condition Number (Can be very expensive !)
+ const bool verbose = true;
+ const unsigned int arppp_a_verbosity_level = 2;
+ const unsigned int pia_verbosity_level = 1;
+ MatrixInfo<MatrixCT> matrixInfo(stiffnessMatrix_,verbose,arppp_a_verbosity_level,pia_verbosity_level);
+ std::cout << "Get condition number of Stiffness_CE: " << matrixInfo.getCond2(true) << std::endl;
+
+ ///////////////////////////////////
+ // --- Choose Solver ---
+ // 1 : CG-Solver
+ // 2 : GMRES
+ // 3 : QR (default)
+ // 4 : UMFPACK
+ ///////////////////////////////////
+ unsigned int Solvertype = parameterSet_.get<unsigned int>("Solvertype", 3);
+ unsigned int Solver_verbosity = parameterSet_.get<unsigned int>("Solver_verbosity", 2);
+
+ // --- set initial values for solver
+ x_1_ = load_alpha1_;
+ x_2_ = load_alpha2_;
+ x_3_ = load_alpha3_;
+
+ Dune::Timer SolverTimer;
+ if (Solvertype==1) // CG - SOLVER
+ {
+ std::cout << "------------ CG - Solver ------------" << std::endl;
+ MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);
+
+ // Sequential incomplete LU decomposition as the preconditioner
+ SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_,1.0);
+ int iter = parameterSet_.get<double>("iterations_CG", 999);
+ // Preconditioned conjugate-gradient solver
+ CGSolver<VectorCT> solver(op,
+ ilu0, //NULL,
+ 1e-8, // desired residual reduction factorlack
+ iter, // maximum number of iterations
+ Solver_verbosity,
+ true // Try to estimate condition_number
+ ); // verbosity of the solver
+ InverseOperatorResult statistics;
+ std::cout << "solve linear system for x_1.\n";
+ solver.apply(x_1_, load_alpha1_, statistics);
+ std::cout << "solve linear system for x_2.\n";
+ solver.apply(x_2_, load_alpha2_, statistics);
+ std::cout << "solve linear system for x_3.\n";
+ solver.apply(x_3_, load_alpha3_, statistics);
+ log_ << "Solver-type used: " <<" CG-Solver" << std::endl;
+
+ std::cout << "statistics.converged " << statistics.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statistics.iterations << std::endl;
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+ else if (Solvertype==2) // GMRES - SOLVER
+ {
+ std::cout << "------------ GMRES - Solver ------------" << std::endl;
+ // Turn the matrix into a linear operator
+ MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_);
+
+ // Fancy (but only) way to not have a preconditioner at all
+ Richardson<VectorCT,VectorCT> preconditioner(1.0);
+
+ // Construct the iterative solver
+ RestartedGMResSolver<VectorCT> solver(
+ stiffnessOperator, // Operator to invert
+ preconditioner, // Preconditioner
+ 1e-10, // Desired residual reduction factor
+ 500, // Number of iterations between restarts,
+ // here: no restarting
+ 500, // Maximum number of iterations
+ Solver_verbosity); // Verbosity of the solver
+
+ // Object storing some statistics about the solving process
+ InverseOperatorResult statistics;
+
+ // solve for different Correctors (alpha = 1,2,3)
+ solver.apply(x_1_, load_alpha1_, statistics); //load_alpha1 now contains the corresponding residual!!
+ solver.apply(x_2_, load_alpha2_, statistics);
+ solver.apply(x_3_, load_alpha3_, statistics);
+ log_ << "Solver-type used: " <<" GMRES-Solver" << std::endl;
+
+ std::cout << "statistics.converged " << statistics.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statistics.iterations << std::endl;
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+ else if ( Solvertype==3)// QR - SOLVER
+ {
+ std::cout << "------------ QR - Solver ------------" << std::endl;
+ log_ << "solveLinearSystems: We use QR solver.\n";
+ //TODO install suitesparse
+ SPQR<MatrixCT> sPQR(stiffnessMatrix_);
+ sPQR.setVerbosity(1);
+ InverseOperatorResult statisticsQR;
+
+ sPQR.apply(x_1_, load_alpha1_, statisticsQR);
+ std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+ sPQR.apply(x_2_, load_alpha2_, statisticsQR);
+ std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+ sPQR.apply(x_3_, load_alpha3_, statisticsQR);
+ std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+ log_ << "Solver-type used: " <<" QR-Solver" << std::endl;
+
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+ else if (Solvertype==4)// UMFPACK - SOLVER
+ {
+ std::cout << "------------ UMFPACK - Solver ------------" << std::endl;
+ log_ << "solveLinearSystems: We use UMFPACK solver.\n";
+
+ Dune::Solvers::UMFPackSolver<MatrixCT,VectorCT> solver;
+ solver.setProblem(stiffnessMatrix_,x_1_,load_alpha1_);
+ // solver.preprocess();
+ solver.solve();
+ solver.setProblem(stiffnessMatrix_,x_2_,load_alpha2_);
+ // solver.preprocess();
+ solver.solve();
+ solver.setProblem(stiffnessMatrix_,x_3_,load_alpha3_);
+ // solver.preprocess();
+ solver.solve();
+ // sPQR.apply(x_1, load_alpha1, statisticsQR);
+ // std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ // std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ // std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+ // sPQR.apply(x_2, load_alpha2, statisticsQR);
+ // std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ // std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ // std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
+ // sPQR.apply(x_3, load_alpha3, statisticsQR);
+ // std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
+ // std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
+ // 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;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Extract phi_alpha & M_alpha coefficients
+ ////////////////////////////////////////////////////////////////////////////////////
+ phi_1_.resize(basis_.size());
+ phi_1_ = 0;
+ phi_2_.resize(basis_.size());
+ phi_2_ = 0;
+ phi_3_.resize(basis_.size());
+ phi_3_ = 0;
+
+ for(size_t i=0; i<basis_.size(); i++)
+ {
+ phi_1_[i] = x_1_[i];
+ phi_2_[i] = x_2_[i];
+ phi_3_[i] = x_3_[i];
+ }
+ for(size_t i=0; i<3; i++)
+ {
+ m_1_[i] = x_1_[phiOffset_+i];
+ m_2_[i] = x_2_[phiOffset_+i];
+ m_3_[i] = x_3_[phiOffset_+i];
+ }
+ // assemble M_alpha's (actually iota(M_alpha) )
+
+ // MatrixRT M_1(0), M_2(0), M_3(0);
+
+ M1_ = 0;
+ M2_ = 0;
+ M3_ = 0;
+
+ for(size_t i=0; i<3; i++)
+ {
+ M1_ += m_1_[i]*MatrixBasisContainer_[i];
+ M2_ += m_2_[i]*MatrixBasisContainer_[i];
+ M3_ += m_3_[i]*MatrixBasisContainer_[i];
+ }
+
+ std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
+ printmatrix(std::cout, M1_, "Corrector-Matrix M_1", "--");
+ printmatrix(std::cout, M2_, "Corrector-Matrix M_2", "--");
+ printmatrix(std::cout, M3_, "Corrector-Matrix M_3", "--");
+ log_ << "---------- OUTPUT ----------" << std::endl;
+ log_ << " --------------------" << std::endl;
+ log_ << "Corrector-Matrix M_1: \n" << M1_ << std::endl;
+ log_ << " --------------------" << std::endl;
+ log_ << "Corrector-Matrix M_2: \n" << M2_ << std::endl;
+ log_ << " --------------------" << std::endl;
+ log_ << "Corrector-Matrix M_3: \n" << M3_ << std::endl;
+ log_ << " --------------------" << std::endl;
+
+
+ if(parameterSet_.get<bool>("write_IntegralMean", false))
+ {
+ std::cout << "check integralMean phi_1: " << std::endl;
+ auto A = integralMean(phi_1_);
+ for(size_t i=0; i<3; i++)
+ {
+ std::cout << "Integral-Mean phi_1 : " << A[i] << std::endl;
+ }
+ }
+ if(substract_integralMean)
+ {
+ std::cout << " --- substracting integralMean --- " << std::endl;
+ subtractIntegralMean(phi_1_);
+ subtractIntegralMean(phi_2_);
+ subtractIntegralMean(phi_3_);
+ subtractIntegralMean(x_1_);
+ subtractIntegralMean(x_2_);
+ subtractIntegralMean(x_3_);
+ //////////////////////////////////////////
+ // Check Integral-mean again:
+ //////////////////////////////////////////
+ if(parameterSet_.get<bool>("write_IntegralMean", false))
+ {
+ auto A = integralMean(phi_1_);
+ for(size_t i=0; i<3; i++)
+ {
+ std::cout << "Integral-Mean phi_1 (Check again) : " << A[i] << std::endl;
+ }
+ }
+ }
+ /////////////////////////////////////////////////////////
+ // Write Solution (Corrector Coefficients) in Logs
+ /////////////////////////////////////////////////////////
+ if(parameterSet_.get<bool>("write_corrector_phi1", false))
+ {
+ log_ << "\nSolution of Corrector problems:\n";
+ log_ << "\n Corrector_phi1:\n";
+ log_ << x_1_ << std::endl;
+ }
+ if(parameterSet_.get<bool>("write_corrector_phi2", false))
+ {
+ log_ << "-----------------------------------------------------";
+ log_ << "\n Corrector_phi2:\n";
+ log_ << x_2_ << std::endl;
+ }
+ if(parameterSet_.get<bool>("write_corrector_phi3", false))
+ {
+ log_ << "-----------------------------------------------------";
+ log_ << "\n Corrector_phi3:\n";
+ log_ << x_3_ << std::endl;
+ }
+
+ }
+
+
+ // -----------------------------------------------------------------
+ // --- Write Correctos to VTK:
+ void writeCorrectorsVTK(const int level)
+ {
+ std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/CellProblem-result";
+ std::cout << "vtkOutputName:" << vtkOutputName << std::endl;
+
+ VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
+ vtkWriter.addVertexData(
+ Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_1_),
+ VTK::FieldInfo("Corrector phi_1 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(
+ Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_2_),
+ VTK::FieldInfo("Corrector phi_2 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.addVertexData(
+ Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_3_),
+ VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
+ vtkWriter.write(vtkOutputName + "-level"+ std::to_string(level));
+ std::cout << "wrote Corrector-VTK data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;
+ }
+
+ // -----------------------------------------------------------------
+ // --- Compute norms of the corrector functions:
+ void computeNorms()
+ {
+ computeL2Norm();
+ computeL2SymGrad();
+
+ std::cout<< "Frobenius-Norm of M1_: " << M1_.frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M2_: " << M2_.frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M3_: " << M3_.frobenius_norm() << std::endl;
+ }
+
+ void computeL2Norm()
+ {
+ // IMPLEMENTATION with makeDiscreteGlobalBasisFunction
+ double error_1 = 0.0;
+ double error_2 = 0.0;
+ double error_3 = 0.0;
+
+ auto localView = basis_.localView();
+ auto GVFunc_1 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_);
+ auto GVFunc_2 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_);
+ auto GVFunc_3 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_);
+ auto localfun_1 = localFunction(GVFunc_1);
+ auto localfun_2 = localFunction(GVFunc_2);
+ auto localfun_3 = localFunction(GVFunc_3);
+
+ for(const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ localfun_1.bind(element);
+ localfun_2.bind(element);
+ localfun_3.bind(element);
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+
+ for(const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = element.geometry().integrationElement(quadPos);
+ error_1 += localfun_1(quadPos)*localfun_1(quadPos) * quadPoint.weight() * integrationElement;
+ error_2 += localfun_2(quadPos)*localfun_2(quadPos) * quadPoint.weight() * integrationElement;
+ error_3 += localfun_3(quadPos)*localfun_3(quadPos) * quadPoint.weight() * integrationElement;
+ }
+ }
+ std::cout << "L2-Norm(Corrector phi_1): " << sqrt(error_1) << std::endl;
+ std::cout << "L2-Norm(Corrector phi_2): " << sqrt(error_2) << std::endl;
+ std::cout << "L2-Norm(Corrector phi_3): " << sqrt(error_3) << std::endl;
+
+ return;
+ }
+
+ void computeL2SymGrad()
+ {
+ double error_1 = 0.0;
+ double error_2 = 0.0;
+ double error_3 = 0.0;
+
+ auto localView = basis_.localView();
+
+ auto GVFunc_1 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
+ auto GVFunc_2 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
+ auto GVFunc_3 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
+ auto localfun_1 = localFunction(GVFunc_1);
+ auto localfun_2 = localFunction(GVFunc_2);
+ auto localfun_3 = localFunction(GVFunc_3);
+
+ for (const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ localfun_1.bind(element);
+ localfun_2.bind(element);
+ localfun_3.bind(element);
+ auto geometry = element.geometry();
+
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
+
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+ for (const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ const auto integrationElement = geometry.integrationElement(quadPos);
+
+ auto scaledSymGrad1 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_1(quadPos)));
+ auto scaledSymGrad2 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_2(quadPos)));
+ auto scaledSymGrad3 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_3(quadPos)));
+
+ error_1 += scalarProduct(scaledSymGrad1,scaledSymGrad1) * quadPoint.weight() * integrationElement;
+ error_2 += scalarProduct(scaledSymGrad2,scaledSymGrad2) * quadPoint.weight() * integrationElement;
+ error_3 += scalarProduct(scaledSymGrad3,scaledSymGrad3) * quadPoint.weight() * integrationElement;
+ }
+ }
+ std::cout << "L2-Norm(Symmetric scaled gradient of Corrector phi_1): " << sqrt(error_1) << std::endl;
+ std::cout << "L2-Norm(Symmetric scaled gradient of Corrector phi_2): " << sqrt(error_2) << std::endl;
+ std::cout << "L2-Norm(Symmetric scaled gradient of Corrector phi_3): " << sqrt(error_3) << std::endl;
+ return;
+ }
+
+
+
+ // -----------------------------------------------------------------
+ // --- Check Orthogonality relation Paper (75)
+ auto check_Orthogonality()
+ {
+ std::cout << "Check Orthogonality..." << std::endl;
+
+ auto localView = basis_.localView();
+
+ auto GVFunc_1 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
+ auto GVFunc_2 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
+ auto GVFunc_3 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
+ auto localfun_1 = localFunction(GVFunc_1);
+ auto localfun_2 = localFunction(GVFunc_2);
+ auto localfun_3 = localFunction(GVFunc_3);
+ const std::array<decltype(localfun_1)*,3> phiDerContainer = {&localfun_1 , &localfun_2 , &localfun_3 };
+
+ auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+
+
+
+ auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ auto mu = localFunction(muGridF);
+ auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ auto lambda= localFunction(lambdaGridF);
+
+ for(int a=0; a<3; a++)
+ for(int b=0; b<3; b++)
+ {
+ double energy = 0.0;
+
+ auto DerPhi1 = *phiDerContainer[a];
+ auto DerPhi2 = *phiDerContainer[b];
+
+ auto matrixFieldGGVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer_[a], basis_.gridView());
+ auto matrixFieldG = localFunction(matrixFieldGGVF);
+ // auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
+ // auto matrixFieldB = localFunction(matrixFieldBGVF);
+
+ // using GridView = typename Basis::GridView;
+ // using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ // using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+ // std::cout << "Press key.." << std::endl;
+ // std::cin.get();
+
+ // TEST
+ // FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
+ // printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
+
+ for (const auto& e : elements(basis_.gridView()))
+ {
+ localView.bind(e);
+ matrixFieldG.bind(e);
+ DerPhi1.bind(e);
+ DerPhi2.bind(e);
+ mu.bind(e);
+ lambda.bind(e);
+ localIndicatorFunction.bind(e);
+
+ double elementEnergy = 0.0;
+ //double elementEnergy_HP = 0.0;
+
+ auto geometry = e.geometry();
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ // int orderQR = 0;
+ // int orderQR = 1;
+ // int orderQR = 2;
+ // int orderQR = 3;
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
+
+ for (const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = geometry.integrationElement(quadPos);
+
+
+ auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + *mContainer[b];
+
+ auto strain1 = DerPhi1(quadPos);
+
+
+ strain1 = crossSectionDirectionScaling(1.0/gamma_, strain1);
+ strain1 = sym(strain1);
+
+
+ auto G = matrixFieldG(quadPos);
+ // auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
+
+ auto tmp = G + *mContainer[a] + strain1;
+
+ double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
+
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+ // elementEnergy += strain1 * quadPoint.weight() * integrationElement;
+ //elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
+ }
+ energy += elementEnergy;
+ //higherPrecEnergy += elementEnergy_HP;
+ }
+ if(parameterSet_.get<bool>("print_debug", false))
+ std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl;
+
+ if(energy > 1e-1)
+ std::cout << "WARNING: check_Orthogonality failed! apparently orthogonality (75) not satisfied on discrete level" << std::endl;
+
+ }
+ return 0;
+ }
+
+
+
+ // --- Check wheter stiffness matrix is symmetric
+ void checkSymmetry()
+ {
+ std::cout << " Check Symmetry of stiffness matrix..." << std::endl;
+
+ auto localView = basis_.localView();
+
+ for (const auto& element : elements(basis_.gridView()))
+ {
+ localView.bind(element);
+ const int localPhiOffset = localView.size();
+
+ for(size_t i=0; i<localPhiOffset; i++)
+ for(size_t j=0; j<localPhiOffset; j++ )
+ {
+ auto row = localView.index(i);
+ auto col = localView.index(j);
+ if(abs( stiffnessMatrix_[row][col] - stiffnessMatrix_[col][row]) > 1e-12 )
+ std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
+ }
+ for(size_t i=0; i<localPhiOffset; i++)
+ for(size_t m=0; m<3; m++)
+ {
+ auto row = localView.index(i);
+ if(abs( stiffnessMatrix_[row][phiOffset_+m] - stiffnessMatrix_[phiOffset_+m][row]) > 1e-12 )
+ std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
+
+ }
+ for(size_t m=0; m<3; m++ )
+ for(size_t n=0; n<3; n++ )
+ {
+ if(abs(stiffnessMatrix_[phiOffset_+m][phiOffset_+n] - stiffnessMatrix_[phiOffset_+n][phiOffset_+m]) > 1e-12 )
+ std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
+ }
+ }
+ std::cout << "--- Symmetry test passed ---" << std::endl;
+ }
+
+
+
+
+
+}; // end class
+
+#endif
\ No newline at end of file
diff --git a/src/deprecated_code/30-8-22_TestingFiles-Comments/prestrainedMaterial.hh b/src/deprecated_code/30-8-22_TestingFiles-Comments/prestrainedMaterial.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3a274edcf111027ab41673537d10aeea589c6d9b
--- /dev/null
+++ b/src/deprecated_code/30-8-22_TestingFiles-Comments/prestrainedMaterial.hh
@@ -0,0 +1,284 @@
+#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
+
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/functions/gridfunctions/gridviewentityset.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/fufem/dunepython.hh>
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/materialDefinitions.hh>
+
+
+using namespace Dune;
+using namespace MatrixOperations;
+using std::pow;
+using std::abs;
+using std::sqrt;
+using std::sin;
+using std::cos;
+
+using std::shared_ptr;
+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();
+// }
+
+template <class GridView> // needed for GridViewFunctions
+class prestrainedMaterial
+{
+
+public:
+ static const int dimworld = 3; //GridView::dimensionworld;
+ static const int dim = 3; //const int dim = Domain::dimension;
+
+ 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>;
+ using FuncScalar = std::function< double(const Domain&) >;
+ using Func2int = std::function< int(const Domain&) >;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
+ using Func2TensorPhase = std::function< MatrixRT(const MatrixRT& ,const int&) >;
+ using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
+ using MatrixPhase = std::function< MatrixRT(const int&) >;
+
+protected:
+
+ const GridView& gridView_;
+ const ParameterTree& parameterSet_;
+ std::string materialFunctionName_;
+
+ // MatrixFunc L1_;
+ // MatrixFunc L2_;
+ // MatrixFunc L3_;
+ // const FuncScalar mu_;
+ // const FuncScalar lambda_;
+ // const int phases_;
+ // Func2Tensor elasticityTensor_;
+ // Func2TensorParam elasticityTensor_;
+ Func2TensorPhase elasticityTensor_;
+ MatrixPhase prestrain_;
+
+
+ Func2int indicatorFunction_;
+ 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)
+ : gridView_(gridView),
+ parameterSet_(parameterSet)
+ {
+ std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material_1");
+ std::cout << "materialFunctionName_ : " << materialFunctionName_ << std::endl;
+
+ setupMaterial(materialFunctionName_);
+
+ indicatorFunctionGVF_ = Dune::Functions::makeGridViewFunction(indicatorFunction_, gridView_);
+ localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
+
+ // 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_);
+ // L1_ = Python::make_function<MatrixRT>(module.get("L1"));
+ // L2_ = Python::make_function<MatrixRT>(module.get("L2"));
+ // L3_ = Python::make_function<MatrixRT>(module.get("L3"));
+ // bool isotropic_ = true; // read from module File
+ // auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
+ }
+
+ //---function that determines elasticity Tensor
+ void setupMaterial(const std::string name) // cant use materialFunctionName_ here!?
+ {
+ // std::cout << "Using material definition:" << name << std::endl;
+ if(name == "material_neukamm")
+ {
+ indicatorFunction_ = indicatorFunction_material_neukamm<Domain>;
+ elasticityTensor_ = material_neukamm;
+ prestrain_ = prestrain_material_neukamm;
+ }
+ else if(name == "two_phase_material_1")
+ {
+ indicatorFunction_ = indicatorFunction_two_phase_material_1<Domain>;
+ elasticityTensor_ = two_phase_material_1;
+ prestrain_ = prestrain_two_phase_material_1;
+ }
+ else if(name == "two_phase_material_2")
+ {
+ indicatorFunction_ = indicatorFunction_two_phase_material_2<Domain>;
+ elasticityTensor_ = two_phase_material_2;
+ prestrain_ = prestrain_two_phase_material_2;
+ }
+ else if(name == "homogeneous")
+ {
+ indicatorFunction_ = indicatorFunction_material_homogeneous<Domain>;
+ elasticityTensor_ = material_homogeneous;
+ prestrain_ = prestrain_homogeneous;
+ }
+ else
+ DUNE_THROW(Exception, "There exists no material with this name ");
+ }
+
+
+ //--- apply elasticityTensor_ to input Matrix G at position x
+ // input: Matrix G, material-phase
+ MatrixRT ElasticityTensor(const MatrixRT& G, const int& phase) const
+ {
+ return elasticityTensor_(G,phase);
+
+ }
+
+ int IndicatorFunction(const Domain& x) const
+ {
+ return indicatorFunction_(x);
+ }
+
+ // 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;
+ // std::cout << "-1/2" << -1/2 << std::endl;
+ // if ((x[1]< (-0.5+theta)) and (x[2]>(0.5-theta)))
+ // {
+ // return 2; //#Phase1
+ // }
+ // else if ((x[0] < (-0.5+theta)) and (x[2]<(-0.5+theta)))
+ // {
+ // return 1; //#Phase2
+ // }
+ // else
+ // {
+ // return 0; //#Phase3
+ // }
+ // }
+
+
+
+ // 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));
+ // }
+
+
+
+
+ // 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 //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_);
+ // // auto IFunction = this->getIndicatorFunction();
+ // // 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);
+ // // auto tmp = elasticityTensor_(G,phase);
+ // // auto tmp = material_1(G,indicatorFunction_(x));
+ // // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
+ // // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
+ // // return tmp;
+ // // return material_1(G,indicatorFunction_(x));
+ // }
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+ // MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
+ // {
+ // //--- apply elasticityTensor_ to input Matrix G at position x
+ // return elasticityTensor_(G,x);
+
+ // }
+
+ // MatrixRT applyElasticityTensorLocal(const MatrixRT& G, const Domain& x) const
+ // {
+ // //--- apply elasticityTensor_ to input Matrix G at position x (local coordinates)
+ // // MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
+
+ // if (indicatorFunction_(x) == 1)
+ // return L1_(G);
+ // else if (indicatorFunction_(x) == 2)
+ // return L2_(G);
+ // else
+ // return L3_(G);
+
+ // }
+
+
+ // -----------------------------------------------------------------
+ // --- write material (grid)functions to VTK
+ void write_materialFunctions()
+ {
+
+
+ return;
+
+ };
+ ///////////////////////////////
+ // Getter
+ ///////////////////////////////
+ ParameterTree getParameterSet() const {return parameterSet_;}
+ // Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
+ // Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
+ Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
+ MatrixPhase getPrestrainFunction() const {return prestrain_;}
+
+ 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_);}
+
+ //getLameParameters() .. Throw Exception if isotropic_ = false!
+
+
+
+
+}; // end class
+
+
+#endif
diff --git a/src/Cell-Problem_muGamma.cc b/src/deprecated_code/Cell-Problem_muGamma.cc
similarity index 100%
rename from src/Cell-Problem_muGamma.cc
rename to src/deprecated_code/Cell-Problem_muGamma.cc