diff --git a/dune/microstructure/prestrain_material_geometry.hh b/dune/microstructure/prestrain_material_geometry.hh
index 092fa5474c74e9cc72e6ab4ad89121e77c67dbbf..ae4560d6e10c7a366d09b772eee5324544f1e208 100644
--- a/dune/microstructure/prestrain_material_geometry.hh
+++ b/dune/microstructure/prestrain_material_geometry.hh
@@ -1,5 +1,5 @@
-#ifndef DUNE_MICROSTRUCTURE_PRESTRAI_MATERIAL_GEOMETRY_HH
-#define DUNE_MICROSTRUCTURE_PRESTRAI_MATERIAL_GEOMETRY_HH
+#ifndef DUNE_MICROSTRUCTURE_PRESTRAIN_MATERIAL_GEOMETRY_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAIN_MATERIAL_GEOMETRY_HH
#include <dune/grid/uggrid.hh>
@@ -18,16 +18,14 @@ using std::cos;
-
+template <int dim>
class IsotropicMaterialImp
{
public:
- static const int dim = 3;
- static const int dimWorld = 3;
+ static const int dimWorld = dim;
using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- using GridView = CellGridType::LeafGridView;
- using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using FuncScalar = std::function< double(const Domain&) >;
@@ -757,16 +755,17 @@ public:
// TODO add log here?
-
+template <int dim>
class PrestrainImp
{
public:
- static const int dim = 3;
+// static const int dim = 3;
static const int dimWorld = 3;
using CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
- using GridView = CellGridType::LeafGridView;
- using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+// using GridView = CellGridType::LeafGridView;
+ using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+// using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
using ScalarRT = double;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
diff --git a/dune/microstructure/vtk_filler.hh b/dune/microstructure/vtk_filler.hh
index 0443c078b83ec1cf155b6e6e0063874f29a6be46..aa61fb2c1c8b475338d0d9def75e57c8709180ca 100644
--- a/dune/microstructure/vtk_filler.hh
+++ b/dune/microstructure/vtk_filler.hh
@@ -3,7 +3,7 @@
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
-//#include <dune/functions/gridfunctions/gridviewfunction.hh>
+//#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/microstructure/prestrain_material_geometry.hh>
@@ -21,7 +21,7 @@
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> >;
@@ -30,8 +30,8 @@
*/
/*
template<class Basis>
- static auto scalarDiscGlobalBasisFunc(
- const Basis& feBasis,
+ static auto scalarDiscGlobalBasisFunc(
+ const Basis& feBasis,
std::function< double(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& scalarFunc)
{
ScalarCT values;
@@ -43,7 +43,7 @@
template<class Basis>
static auto vectorDiscGlobalBasisFunc(
- const Basis& feBasis,
+ const Basis& feBasis,
std::function< FieldVector< double, dimworld>(const typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate&) >& vectorFunc)
{
VectorCT values;
@@ -76,17 +76,17 @@ public:
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);
-
+ FuncScalar normB = PrestrainImp<dim>().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);
@@ -94,13 +94,13 @@ public:
(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(
@@ -118,8 +118,8 @@ public:
}
- void vtkMaterialCell(const GV& gridView_CE,
- const FuncScalar& mu_Func,
+ void vtkMaterialCell(const GV& gridView_CE,
+ const FuncScalar& mu_Func,
std::string filename)
{
@@ -135,15 +135,15 @@ public:
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,
+
+ void vtkMaterial3DRod(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncScalar& mu_Func,
std::string filename)
{
@@ -168,15 +168,15 @@ public:
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,
+ void vtkProblemCell(const GV& gridView_CE,
+ const FuncMatrix& B_Func, const FuncScalar& mu_Func,
std::string filename)
{
@@ -187,7 +187,7 @@ public:
// 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];
@@ -218,10 +218,10 @@ public:
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);
+ FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
ScalarCT prestrain_Coeff;
Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
@@ -233,14 +233,14 @@ public:
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,
+ void vtkProblemCellBVec(const GV& gridView_CE,
+ const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
std::string filename)
{
@@ -251,7 +251,7 @@ public:
// 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];
@@ -289,12 +289,12 @@ public:
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);
+ FuncScalar normB_Func = PrestrainImp<dim>().getNormPrestrain(B_Func);
ScalarCT prestrain_Coeff;
Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
@@ -306,15 +306,15 @@ public:
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,
+ void vtkProblem3DRod(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
std::string filename)
{
@@ -331,7 +331,7 @@ public:
(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]};
@@ -368,11 +368,11 @@ public:
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_Func = PrestrainImp<dim>().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]};
@@ -389,14 +389,14 @@ public:
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,
+ 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)
{
@@ -413,7 +413,7 @@ public:
(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]};
@@ -460,11 +460,11 @@ public:
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_Func = PrestrainImp<dim>().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]};
@@ -486,7 +486,7 @@ public:
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";
@@ -496,9 +496,9 @@ public:
/*
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,
+ 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>;
@@ -514,7 +514,7 @@ public:
// 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];
@@ -547,11 +547,11 @@ public:
(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);
+ FuncScalar normB_Func = getNormPrestrain(B_Func);
ScalarCT prestrain_Coeff;
Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
@@ -563,14 +563,14 @@ public:
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,
+ void vtkProblemCellBVec(const GV& gridView_CE,
+ const FuncMatrix& B_Func, const FuncVector& B_vec_Func, const FuncScalar& mu_Func,
std::string filename)
{
@@ -581,7 +581,7 @@ public:
// 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];
@@ -619,12 +619,12 @@ public:
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);
+ FuncScalar normB_Func = getNormPrestrain(B_Func);
ScalarCT prestrain_Coeff;
Functions::interpolate(scalarFeBasis, prestrain_Coeff, normB_Func);
auto normB_DGBF = Functions::makeDiscreteGlobalBasisFunction< ScalarRT >
@@ -636,15 +636,15 @@ public:
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,
+ void vtkProblem3DRod(const GV& gridView_R3D, double eps,
+ const FuncVector& domain_Func, const FuncMatrix& B_Func, const FuncScalar& mu_Func,
std::string filename)
{
@@ -661,7 +661,7 @@ public:
(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]};
@@ -698,11 +698,11 @@ public:
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_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]};
@@ -719,14 +719,14 @@ public:
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,
+ 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)
{
@@ -743,7 +743,7 @@ public:
(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]};
@@ -790,11 +790,11 @@ public:
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_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]};
@@ -816,7 +816,7 @@ public:
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";
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index dbc2cfd525c614b7137e331113e986d0742a8796..bf2c1f8b7b83b1f1053cbfd966cbc925b08ee502 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -168,12 +168,16 @@ void computeElementStiffnessMatrix(const LocalView& localView,
//////////////////////////////////////////////
MatrixRT G_1 {{1.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}};
- MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.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);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+
+
+// std::cout << "Print QuadratureOrder:" << orderQR << std::endl; //TEST
for (const auto& quadPoint : quad)
{
@@ -342,7 +346,7 @@ void computeElementLoadVector( const LocalView& localView,
//////////////////////////////////////////////
MatrixRT G_1 {{1.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}};
- MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.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
@@ -541,7 +545,7 @@ auto energy(const Basis& basis,
auto geometry = e.geometry();
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 + 5 ); // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
@@ -663,7 +667,8 @@ auto integralMean(const Basis& basis,
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);
@@ -747,6 +752,7 @@ double computeL2SymError(const Basis& basis,
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; //TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -818,6 +824,8 @@ double computeL2Norm(const Basis& basis,
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);
@@ -870,13 +878,13 @@ int main(int argc, char *argv[])
///////////////////////////////////
// Get Parameters/Data
///////////////////////////////////
- double gamma = parameterSet.get<double>("gamma",3.0); // ratio dimension reduction to homogenization
+ 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_implementation", "analytical_Example");
+ std::string imp = parameterSet.get<std::string>("material_prestrain_imp", "analytical_Example");
double beta = parameterSet.get<double>("beta",2.0);
double mu1 = parameterSet.get<double>("mu1",10.0);;
double mu2 = beta*mu1;
@@ -901,7 +909,7 @@ int main(int argc, char *argv[])
- auto prestrainImp = PrestrainImp(); //NEW
+ auto prestrainImp = PrestrainImp<dim>(); //NEW
auto B_Term = prestrainImp.getPrestrain(parameterSet);
log << "----- Input Parameters -----: " << std::endl;
@@ -974,7 +982,7 @@ int main(int argc, char *argv[])
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
- auto materialImp = IsotropicMaterialImp();
+ auto materialImp = IsotropicMaterialImp<dim>();
auto muTerm = materialImp.getMu(parameterSet);
auto lambdaTerm = materialImp.getLambda(parameterSet);
@@ -1037,7 +1045,7 @@ int main(int argc, char *argv[])
// First order periodic Lagrange-Basis
auto Basis_CE = makeBasis(
gridView_CE,
- power<dim>(
+ power<dim>( // eig dimworld?!?!
Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
flatLexicographic()));
@@ -1062,7 +1070,7 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
// Define "rhs"-functions
/////////////////////////////////////////////////////////
- Func2Tensor x3G_1 = [] (const Domain& x, double sign=1.0) {
+ 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?
};
@@ -1071,7 +1079,7 @@ int main(int argc, char *argv[])
};
Func2Tensor x3G_3 = [] (const Domain& x) {
- return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ 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}};
};
Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {return -1.0*x3G_1(x);};
@@ -1099,7 +1107,7 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////
MatrixRT G_1 {{1.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}};
- MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.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};
// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
@@ -1277,7 +1285,7 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////
// Substract Integral-mean
- //////////////////////////////////////////////////////////////////////////// // ERROR
+ ////////////////////////////////////////////////////////////////////////////
using SolutionRange = FieldVector<double,dim>;
if(write_IntegralMean)
@@ -1430,7 +1438,9 @@ int main(int argc, char *argv[])
std::cout<< "q1 : " << q1 << std::endl;
std::cout<< "q2 : " << q2 << std::endl;
- std::cout<< "q3 : " << q3 << std::endl;
+ std::cout.precision(10);
+ std::cout << "q3 : " << std::fixed << q3 << std::endl;
+// std::cout<< "q3 : " << q3 << std::endl;
printvector(std::cout, B_hat, "B_hat", "--");
printvector(std::cout, Beff, "Beff", "--");
diff --git a/src/Compute_MuGamma.cc b/src/Compute_MuGamma.cc
index 8482b7dbf91846cf29e968a0b42a06c72dcde7c6..ad14fab74274f0ad10e19b4be13010f1f2eaf54e 100644
--- a/src/Compute_MuGamma.cc
+++ b/src/Compute_MuGamma.cc
@@ -64,19 +64,6 @@ using namespace Dune;
-
-
-
-
-
-
-
-
-
-
-
-
-
template<class LocalView, class Matrix, class LocalFunction>
void assembleElementStiffnessMatrix(const LocalView& localView,
Matrix& elementMatrix,
@@ -702,10 +689,6 @@ double computeMuGamma(const Basis& basis,
*/
-
-
-
-
// Check whether two points are equal on R/Z x R/Z
auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& y)
{
@@ -714,19 +697,6 @@ double computeMuGamma(const Basis& basis,
};
-// // Check whether two points are equal on R/Z x R/Z x R
-// auto equivalent3D = [](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]))
-// );
-// };
-//
-//
-
-
-
int main(int argc, char *argv[])
{
@@ -760,21 +730,14 @@ int main(int argc, char *argv[])
// std::array<int, dim> nElements = {16,16};
// std::array<int, dim> nElements = {8,8};
- /*
- 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, dim> nElements = {16,16,16};*/
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
CellGridType grid_CE(lower,upper,nElements);
using GridView = CellGridType::LeafGridView;
const GridView gridView = grid_CE.leafGridView();
-
-
using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
-
-
+ // ----------- INPUT PARAMETERS -----------------------------
double gamma = parameterSet.get<double>("gamma", 1.0);
double theta = parameterSet.get<double>("theta", 1.0/4.0);
double mu1 = parameterSet.get<double>("mu1", 10.0);
@@ -797,10 +760,12 @@ int main(int argc, char *argv[])
};
-// auto materialImp = IsotropicMaterialImp();
-// auto muTermT = materialImp.getMu(parameterSet);
+ auto materialImp = IsotropicMaterialImp<dim>();
+ auto muTermT = materialImp.getMu(parameterSet);
- auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView);
+// auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView);
+// auto muLocal = localFunction(muGridF);
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTermT, gridView);
auto muLocal = localFunction(muGridF);