#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);
// }
///////////////////////////////
// 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!
///////////////////////////////
// VTK - Writer
///////////////////////////////
// -----------------------------------------------------------------
// --- write material (grid)functions to VTK
void writeVTKMaterialFunctions(std::array<int, dim> nElements, const int level)
{
std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
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},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> indicatorFunction_CoeffP0;
Functions::interpolate(scalarP0FeBasis, indicatorFunction_CoeffP0, indicatorFunction_);
auto indicatorFunction_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, indicatorFunction_CoeffP0);
std::vector<double> indicatorFunction_CoeffP1;
Functions::interpolate(scalarP1FeBasis, indicatorFunction_CoeffP1, indicatorFunction_);
auto indicatorFunction_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, indicatorFunction_CoeffP1);
VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
MaterialVtkWriter.addVertexData(
indicatorFunction_P0,
VTK::FieldInfo("indicatorFunction_P0", VTK::FieldInfo::Type::scalar, 1));
MaterialVtkWriter.addCellData(
indicatorFunction_P1,
VTK::FieldInfo("indicatorFunction_P1", 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;
return;
};
// -----------------------------------------------------------------
// --- write prestrain (grid)functions to VTK
// void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
// {
// std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
// 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},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> prestrainFunction_CoeffP0;
// Functions::interpolate(scalarP0FeBasis, prestrainFunction_CoeffP0, prestrain_);
// auto prestrainFunction_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, prestrainFunction_CoeffP0);
// std::vector<double> prestrainFunction_CoeffP1;
// Functions::interpolate(scalarP1FeBasis, prestrainFunction_CoeffP1, prestrain_);
// auto prestrainFunction_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, prestrainFunction_CoeffP1);
// VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
// MaterialVtkWriter.addVertexData(
// prestrainFunction_P0,
// VTK::FieldInfo("prestrainFunction_P0", VTK::FieldInfo::Type::scalar, 1));
// MaterialVtkWriter.addCellData(
// prestrainFunction_P1,
// VTK::FieldInfo("prestrainFunction_P1", VTK::FieldInfo::Type::scalar, 1));
// MaterialVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
// std::cout << "wrote data to file:" + outputPath + "/PrestrainFunctions-level" + std::to_string(level) << std::endl;
// return;
// };
// void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
// {
// 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;
// // };
}; // end class
#endif