diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index ad03b2ae60670aaca702f35a2c16ec1c9d8af898..2ceacb0e89613a710395c1f3d85031e57f56f682 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -722,7 +722,11 @@ public:
using GridView = CellGridType::LeafGridView;
using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+ using ScalarRT = double;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ using FuncScalar = std::function< ScalarRT(const Domain&) >;
+ using FuncMatrix = std::function< MatrixRT(const Domain&) >;
+ using FuncVector = std::function< VectorRT(const Domain&) >;
protected:
// double p1, p2, theta;
// double width; //Cell geometry
@@ -919,6 +923,17 @@ public:
}
// TODO ANISOTROPIC PRESSURE
}
+
+
+
+ FuncScalar getNormPrestrain(const FuncMatrix& B)
+ {
+ FuncScalar normB = [&](const Domain& z) {
+ return norm(B(z));
+ };
+ return normB;
+ }
+
};
diff --git a/dune/microstructure/vtk_filler.hh b/dune/microstructure/vtk_filler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0443c078b83ec1cf155b6e6e0063874f29a6be46
--- /dev/null
+++ b/dune/microstructure/vtk_filler.hh
@@ -0,0 +1,826 @@
+#ifndef DUNE_MICROSTRUCTURE_VTK_FILLER_HH
+#define DUNE_MICROSTRUCTURE_VTK_FILLER_HH
+
+#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
+#include <dune/functions/functionspacebases/interpolate.hh>
+//#include <dune/functions/gridfunctions/gridviewfunction.hh>
+#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
+#include <dune/microstructure/prestrain_material_geometry.hh>
+
+
+ using std::string;
+ using namespace Dune;
+ using namespace Functions;
+ using namespace Functions::BasisFactory;
+
+ //static const int dimworld = 3;
+
+ //using GT = UGGrid<dim>; //typename MicrostructuredRodGrid::GT<dim>;
+ //using GV = typename GT::LeafGridView;
+ /*
+ using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+ using VectorRT = FieldVector< double, dimworld>;
+ using ScalarRT = double;
+
+
+ using ScalarCT = std::vector<double>;
+ using VectorCT = BlockVector<FieldVector<double,1> >;
+
+ using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
+*/
+ /*
+ template<class Basis>
+ static auto scalarDiscGlobalBasisFunc(
+ const Basis& feBasis,
+ std::function< double(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& scalarFunc)
+ {
+ ScalarCT values;
+ Functions::interpolate(feBasis, values, scalarFunc);
+ auto discGlobalBasisFunc = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (feBasis, values);
+ return discGlobalBasisFunc;
+ }
+
+ template<class Basis>
+ static auto vectorDiscGlobalBasisFunc(
+ const Basis& feBasis,
+ std::function< FieldVector< double, dimworld>(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& vectorFunc)
+ {
+ VectorCT values;
+ Functions::interpolate(feBasis, values, vectorFunc);
+ auto discGlobalBasisFunc = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(feBasis), HierarchicVectorView(values));
+ return discGlobalBasisFunc;
+ }
+*/
+
+template <int dim>
+class FTKfillerContainer {
+
+public:
+
+ static const int dimworld = 3;
+
+// using GT = UGGrid<dim>; //typename MicrostructuredRodGrid::GT<dim>;
+ using GT =YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
+ using GV = typename GT::LeafGridView;
+ using Domain = typename GV::template Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+ using VectorRT = FieldVector< double, dimworld>;
+ using ScalarRT = double;
+ using FuncMatrix = std::function< MatrixRT(const Domain&) >;
+ using FuncVector = std::function< VectorRT(const Domain&) >;
+ using FuncScalar = std::function< ScalarRT(const Domain&) >;
+
+ using ScalarCT = std::vector<double>;
+ using VectorCT = BlockVector<FieldVector<double,1> >;
+ using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
+
+
+ void vtkPrestrainNorm(const GV& gridView, const FuncMatrix& B, string filename) // Updated
+ {
+ VTKWriter<typename GT::LeafGridView> vtkWriter(gridView);
+ Functions::LagrangeBasis<typename GT::LeafGridView, 1> feBasis(gridView);
+ FuncScalar normB = PrestrainImp().getNormPrestrain(B);
+
+// auto discGlobalBasisFunc = scalarDiscGlobalBasisFunc(feBasis, normB);
+ auto discGlobalBasisFunc = Dune::Functions::makeGridViewFunction(normB, gridView);
+// auto muLocal = localFunction(muGridF);
+
+
+ ScalarCT values;
+ Functions::interpolate(feBasis, values, normB);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (feBasis, values);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+ vtkWriter.write(filename);
+
+ }
+
+ void vtkPrestrainDirector(const typename GT::LeafGridView& gridView, const FuncVector& b, string filename)
+ {
+ VTKWriter<typename GT::LeafGridView> vtkWriter(gridView);
+
+ using namespace Functions::BasisFactory;
+
+ auto basis_CE = makeBasis(
+ gridView,
+ power<dimworld>(
+ lagrange<1>(),
+ flatLexicographic()));//flatInterleaved()
+
+ VectorCT values;
+ Functions::interpolate(basis_CE, values, b);
+ auto b_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, values);
+
+ vtkWriter.addVertexData(b_DGBF, VTK::FieldInfo("prestrain_director", VTK::FieldInfo::Type::vector, dimworld));
+ vtkWriter.write(filename);
+ }
+
+
+ void vtkMaterialCell(const GV& gridView_CE,
+ const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkMaterialCell ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_CE);
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addCellData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkMaterialCell done.\n";
+ }
+
+
+ void vtkMaterial3DRod(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblem ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_R3D);
+
+ // Vector data
+ const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+ VectorCT domain_Coeff;
+ Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+ auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+ vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+
+ void vtkProblemCell(const GV& gridView_CE,
+ const FuncMatrix& B_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblemCell ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_CE);
+
+ // Vector data
+ const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+ FuncVector Be1_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[0];
+ };
+ VectorCT Be1_Coeff;
+ Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+ auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+ vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+ FuncVector Be2_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[1];
+ };
+ VectorCT Be2_Coeff;
+ Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+ auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+ vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[2];
+ };
+ VectorCT Be3_Coeff;
+ Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+ auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+ vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+ FuncScalar normB_Func = PrestrainImp().getNormPrestrain(B_Func);
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+ void vtkProblemCellBVec(const GV& gridView_CE,
+ const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblemCellBVec ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_CE);
+
+ // Vector data
+ const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+ FuncVector Be1_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[0];
+ };
+ VectorCT Be1_Coeff;
+ Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+ auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+ vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be2_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[1];
+ };
+ VectorCT Be2_Coeff;
+ Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+ auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+ vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[2];
+ };
+ VectorCT Be3_Coeff;
+ Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+ auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+ vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+ VectorCT B_vec_Coeff;
+ Functions::interpolate(basis_CE, B_vec_Coeff, B_vec_Func);
+ auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, B_vec_Coeff);
+ vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+ FuncScalar normB_Func = PrestrainImp().getNormPrestrain(B_Func);
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+
+ void vtkProblem3DRod(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblem ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_R3D);
+
+ // Vector data
+ const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+ VectorCT domain_Coeff;
+ Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+ auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+ vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[0];
+ };
+ VectorCT Be1_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+ auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+ vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[1];
+ };
+ VectorCT Be2_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+ auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+ vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[2];
+ };
+ VectorCT Be3_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+ auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+ vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+ FuncScalar normB_Func = PrestrainImp().getNormPrestrain(B_Func);
+ FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return normB_Func(x_Ce);
+ };
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+ void vtkProblem3DRodBVec(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblem ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_R3D);
+
+ // Vector data
+ const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+ VectorCT domain_Coeff;
+ Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+ auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+ vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[0];
+ };
+ VectorCT Be1_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+ auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+ vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[1];
+ };
+ VectorCT Be2_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+ auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+ vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[2];
+ };
+ VectorCT Be3_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+ auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+ vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+ FuncVector B_vec_R3D_Func = [B_vec_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_vec_Func(x_Ce);
+ };
+ VectorCT B_vec_Coeff;
+ Functions::interpolate(basis_R3D, B_vec_Coeff, B_vec_R3D_Func);
+ auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_R3D, B_vec_Coeff);
+ vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+ FuncScalar normB_Func = PrestrainImp().getNormPrestrain(B_Func);
+ FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return normB_Func(x_Ce);
+ };
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ FuncScalar mu_R3D_Func = [mu_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return mu_Func(x_Ce);
+ };
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_R3D_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+};
+
+/*
+ template<int dim>
+ void vtkProblemCell(const typename UGGrid<dim>::LeafGridView& gridView_CE,
+ const std::function< MatrixRT(const typename UGGrid<dim>::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate&)& B_Func,
+ const std::function< ScalarRT(const typename UGGrid<dim>::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate&)& mu_Func,
+ std::string filename)
+ {
+ using GT = UGGrid<dim>;
+ using GV = typename GT::LeafGridView;
+ using Domain = typename GV::template Codim<0>::Geometry::GlobalCoordinate;
+ using FuncMatrix = std::function< MatrixRT(const Domain&) >;
+ using FuncVector = std::function< VectorRT(const Domain&) >;
+ using FuncScalar = std::function< ScalarRT(const Domain&) >;
+ std::cout << "vtkProblemCell ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_CE);
+
+ // Vector data
+ const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+ FuncVector Be1_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[0];
+ };
+ VectorCT Be1_Coeff;
+ Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+ auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+ vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be2_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[1];
+ };
+ VectorCT Be2_Coeff;
+ Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+ auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+ vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[2];
+ };
+ VectorCT Be3_Coeff;
+ Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+ auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+ vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+ FuncScalar normB_Func = getNormPrestrain(B_Func);
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+ void vtkProblemCellBVec(const GV& gridView_CE,
+ const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblemCell ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_CE);
+
+ // Vector data
+ const auto basis_CE = makeBasis(gridView_CE, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+
+ FuncVector Be1_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[0];
+ };
+ VectorCT Be1_Coeff;
+ Functions::interpolate(basis_CE, Be1_Coeff, Be1_Func);
+ auto Be1_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be1_Coeff));
+ vtkWriter.addVertexData(Be1_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be2_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[1];
+ };
+ VectorCT Be2_Coeff;
+ Functions::interpolate(basis_CE, Be2_Coeff, Be2_Func);
+ auto Be2_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be2_Coeff));
+ vtkWriter.addVertexData(Be2_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_Func = [B_Func](const auto& x)
+ {
+ return B_Func(x)[2];
+ };
+ VectorCT Be3_Coeff;
+ Functions::interpolate(basis_CE, Be3_Coeff, Be3_Func);
+ auto Be3_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_CE), HierarchicVectorView(Be3_Coeff));
+ vtkWriter.addVertexData(Be3_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+ VectorCT B_vec_Coeff;
+ Functions::interpolate(basis_CE, B_vec_Coeff, B_vec_Func);
+ auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_CE, B_vec_Coeff);
+ vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_CE);
+
+ FuncScalar normB_Func = getNormPrestrain(B_Func);
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+
+ void vtkProblem3DRod(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblem ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_R3D);
+
+ // Vector data
+ const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+ VectorCT domain_Coeff;
+ Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+ auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+ vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[0];
+ };
+ VectorCT Be1_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+ auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+ vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[1];
+ };
+ VectorCT Be2_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+ auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+ vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[2];
+ };
+ VectorCT Be3_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+ auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+ vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+ FuncScalar normB_Func = getNormPrestrain(B_Func);
+ FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return normB_Func(x_Ce);
+ };
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+
+ void vtkProblem3DRodBVec(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
+ std::string filename)
+ {
+
+ std::cout << "vtkProblem ...\n";
+
+ VTKWriter<GV> vtkWriter(gridView_R3D);
+
+ // Vector data
+ const auto basis_R3D = makeBasis(gridView_R3D, power<dimworld>(lagrange<1>(), flatLexicographic())); //flatInterleaved()
+
+ VectorCT domain_Coeff;
+ Functions::interpolate(basis_R3D, domain_Coeff, domain_Func);
+ auto domain_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(domain_Coeff));
+ vtkWriter.addVertexData(domain_DGBF, VTK::FieldInfo("Domain", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be1_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[0];
+ };
+ VectorCT Be1_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be1_R3D_Coeff, Be1_R3D_Func);
+ auto Be1_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be1_R3D_Coeff));
+ vtkWriter.addVertexData(Be1_R3D_DGBF, VTK::FieldInfo("B_e1", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be2_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[1];
+ };
+ VectorCT Be2_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be2_R3D_Coeff, Be2_R3D_Func);
+ auto Be2_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be2_R3D_Coeff));
+ vtkWriter.addVertexData(Be2_R3D_DGBF, VTK::FieldInfo("B_e2", VTK::FieldInfo::Type::vector, 3));
+
+
+ FuncVector Be3_R3D_Func = [B_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_Func(x_Ce)[2];
+ };
+ VectorCT Be3_R3D_Coeff;
+ Functions::interpolate(basis_R3D, Be3_R3D_Coeff, Be3_R3D_Func);
+ auto Be3_R3D_DGBF = Functions::makeDiscreteGlobalBasisFunction<VectorRT>
+ (Functions::subspaceBasis(basis_R3D), HierarchicVectorView(Be3_R3D_Coeff));
+ vtkWriter.addVertexData(Be3_R3D_DGBF, VTK::FieldInfo("B_e3", VTK::FieldInfo::Type::vector, 3));
+
+ FuncVector B_vec_R3D_Func = [B_vec_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return B_vec_Func(x_Ce);
+ };
+ VectorCT B_vec_Coeff;
+ Functions::interpolate(basis_R3D, B_vec_Coeff, B_vec_R3D_Func);
+ auto B_vec_DGBF = Functions::makeDiscreteGlobalBasisFunction< VectorRT > (basis_R3D, B_vec_Coeff);
+ vtkWriter.addVertexData(B_vec_DGBF, VTK::FieldInfo("B_vec", VTK::FieldInfo::Type::vector, dimworld));
+
+
+
+ // Scalar data
+ Functions::LagrangeBasis<GV, 1> scalarFeBasis(gridView_R3D);
+
+ FuncScalar normB_Func = getNormPrestrain(B_Func);
+ FuncScalar normB_R3D_Func = [normB_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return normB_Func(x_Ce);
+ };
+ ScalarCT prestrain_Coeff;
+ Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_R3D_Func);
+ auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, prestrain_Coeff);
+ vtkWriter.addVertexData(normB_DGBF, VTK::FieldInfo("normB", VTK::FieldInfo::Type::scalar, 1));
+
+ FuncScalar mu_R3D_Func = [mu_Func, eps](const auto& x)
+ {
+ auto x_Ce = {x[0]/eps - int(x[0]/eps),x[1],x[2]};
+ return mu_Func(x_Ce);
+ };
+ ScalarCT material_Coeff;
+ Functions::interpolate(scalarFeBasis, material_Coeff, mu_R3D_Func);
+ auto material_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
+ (scalarFeBasis, material_Coeff);
+ vtkWriter.addVertexData(material_DGBF, VTK::FieldInfo("mu", VTK::FieldInfo::Type::scalar, 1));
+
+ vtkWriter.write(filename);
+
+ std::cout << "vtkProblem done.\n";
+ }
+*/
+
+#endif
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 8b4e0a04611172da523a76d2532678f994d9862a..80d292950162b649d883a1703caa84938441d6e1 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -990,7 +990,11 @@ int main(int argc, char *argv[])
// Print Options
bool print_debug = parameterSet.get<bool>("print_debug", false);
+
bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
+ bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false);
+
+
bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
bool write_corrector_phi2 = parameterSet.get<bool>("write_corrector_phi2", false);
bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
@@ -1550,7 +1554,7 @@ int main(int argc, char *argv[])
- if (write_materialFunctions )
+ if (write_materialFunctions)
{
using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
@@ -1595,10 +1599,24 @@ int main(int argc, char *argv[])
MaterialVtkWriter.write("MaterialFunctions-level"+ std::to_string(level) );
std::cout << "wrote data to file: 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});
+ 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");;
+ }
-
levelCounter++;
} // Level-Loop End