From 78bfc1e8f6a6e2353bb46f872720675096ae5a71 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Sun, 28 Aug 2022 10:35:06 +0200
Subject: [PATCH] Python ElasticityTensor works but is much slower
---
dune/microstructure/CorrectorComputer.hh | 44 +-
.../EffectiveQuantitiesComputer.hh | 33 +-
dune/microstructure/matrix_operations.hh | 21 +
dune/microstructure/prestrainedMaterial.hh | 57 +-
geometries/material.py | 180 ++-
geometries/python_matrix_operations.py | 46 +
src/Cell-Problem-New.cc | 111 +-
.../Cell-Problem-New.cc | 766 ++++++++++
.../CorrectorComputer.hh | 1325 +++++++++++++++++
.../EffectiveQuantitiesComputer.hh | 469 ++++++
.../material.py | 138 ++
.../prestrainedMaterial.hh | 199 +++
12 files changed, 3300 insertions(+), 89 deletions(-)
create mode 100644 geometries/python_matrix_operations.py
create mode 100644 src/deprecated_code/elasticityTensor-globalFunctionVersion/Cell-Problem-New.cc
create mode 100644 src/deprecated_code/elasticityTensor-globalFunctionVersion/CorrectorComputer.hh
create mode 100644 src/deprecated_code/elasticityTensor-globalFunctionVersion/EffectiveQuantitiesComputer.hh
create mode 100644 src/deprecated_code/elasticityTensor-globalFunctionVersion/material.py
create mode 100644 src/deprecated_code/elasticityTensor-globalFunctionVersion/prestrainedMaterial.hh
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 9fdf01aa..17809b83 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -21,7 +21,7 @@ using std::make_shared;
using std::fstream;
-template <class Basis> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
+template <class Basis, class Material> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
class CorrectorComputer {
public:
@@ -49,6 +49,9 @@ protected:
//private:
const Basis& basis_;
+
+ const Material& material_;
+
fstream& log_; // Output-log
const ParameterTree& parameterSet_;
@@ -95,12 +98,14 @@ 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),
@@ -124,7 +129,9 @@ public:
// --- 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_);
@@ -256,6 +263,9 @@ public:
elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
elementMatrix = 0;
+
+ auto elasticityTensor = material_.getElasticityTensor();
+
// LocalBasis-Offset
const int localPhiOffset = localView.size();
@@ -326,7 +336,17 @@ public:
// printmatrix(std::cout, defGradientU , "defGradientU", "--");
defGradientU = crossSectionDirectionScaling((1.0/gamma_),defGradientU);
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+
+ // 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);
+ // 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..
@@ -338,7 +358,10 @@ public:
// "m*phi" & "phi*m" - part
for( size_t m=0; m<3; m++)
{
- double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+ double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),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;
@@ -360,7 +383,11 @@ public:
// std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
// std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
- double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], 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(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)
@@ -466,7 +493,10 @@ public:
defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), 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
@@ -478,7 +508,9 @@ public:
// "f*m"-part
for (size_t m=0; m<3; m++)
{
- double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),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 = 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
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index d4f4b254..6c9d38ce 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -21,9 +21,9 @@ using std::cout;
using std::endl;
// template <class Basis>
-// class EffectiveQuantitiesComputer : public CorrectorComputer<Basis> {
+// class EffectiveQuantitiesComputer : public CorrectorComputer<Basis,Material> {
-template <class Basis>
+template <class Basis, class Material>
class EffectiveQuantitiesComputer {
public:
@@ -32,22 +32,23 @@ public:
static const int dim = Basis::GridView::dimension;
- using Domain = typename CorrectorComputer<Basis>::Domain;
+ using Domain = typename CorrectorComputer<Basis,Material>::Domain;
- using VectorRT = typename CorrectorComputer<Basis>::VectorRT;
- using MatrixRT = typename CorrectorComputer<Basis>::MatrixRT;
+ using VectorRT = typename CorrectorComputer<Basis,Material>::VectorRT;
+ using MatrixRT = typename CorrectorComputer<Basis,Material>::MatrixRT;
- using Func2Tensor = typename CorrectorComputer<Basis>::Func2Tensor;
- using FuncVector = typename CorrectorComputer<Basis>::FuncVector;
+ using Func2Tensor = typename CorrectorComputer<Basis,Material>::Func2Tensor;
+ using FuncVector = typename CorrectorComputer<Basis,Material>::FuncVector;
- using VectorCT = typename CorrectorComputer<Basis>::VectorCT;
+ using VectorCT = typename CorrectorComputer<Basis,Material>::VectorCT;
- using HierarchicVectorView = typename CorrectorComputer<Basis>::HierarchicVectorView;
+ using HierarchicVectorView = typename CorrectorComputer<Basis,Material>::HierarchicVectorView;
protected:
- CorrectorComputer<Basis>& correctorComputer_;
+ CorrectorComputer<Basis,Material>& correctorComputer_;
Func2Tensor prestrain_;
+ const Material& material_;
public:
VectorCT B_load_TorusCV_; //<B, Chi>_L2
@@ -93,10 +94,14 @@ public:
///////////////////////////////
// constructor
///////////////////////////////
- // EffectiveQuantitiesComputer(CorrectorComputer<Basis>& correctorComputer, Func2Tensor prestrain)
+ // EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer, Func2Tensor prestrain)
// : correctorComputer_(correctorComputer), prestrain_(prestrain)
- EffectiveQuantitiesComputer(CorrectorComputer<Basis>& correctorComputer, Func2Tensor prestrain)
- : correctorComputer_(correctorComputer), prestrain_(prestrain)
+ EffectiveQuantitiesComputer(CorrectorComputer<Basis,Material>& correctorComputer,
+ Func2Tensor prestrain,
+ const Material& material)
+ : correctorComputer_(correctorComputer),
+ prestrain_(prestrain),
+ material_(material)
{
// computePrestressLoadCV();
@@ -118,7 +123,7 @@ public:
///////////////////////////////
// getter
///////////////////////////////
- CorrectorComputer<Basis> getCorrectorComputer(){return correctorComputer_;}
+ CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}
const shared_ptr<Basis> getBasis()
{
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index ed403604..33b1e590 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -179,6 +179,27 @@ namespace MatrixOperations {
+
+ extern "C"
+ {
+
+ MatrixRT new_sym(MatrixRT M)
+ {
+ return sym(M);
+ }
+
+
+
+
+
+ }
+
+
+
+
+
+
+
/*
template<double phi>
static bool isInRotatedPlane(double x1, double x2){
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index df3a516d..e255a6b2 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -5,6 +5,9 @@
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/microstructure/matrix_operations.hh>
+#include <dune/functions/gridfunctions/gridviewfunction.hh>
+
+#include <dune/common/parametertree.hh>
#include <dune/fufem/dunepython.hh>
@@ -40,6 +43,7 @@ public:
using FuncScalar = std::function< double(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
+ using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
protected:
@@ -47,6 +51,14 @@ protected:
const GridView& gridView_;
const ParameterTree& parameterSet_;
+ int Phases_;
+
+ MatrixFunc L1_;
+ MatrixFunc L2_;
+ MatrixFunc L3_;
+
+
+
// const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
@@ -64,6 +76,11 @@ protected:
// Func2Tensor elasticityTensor_;
Func2TensorParam elasticityTensor_;
+ // FuncScalar indicatorFunction_;
+
+ GridViewFunction<double(const Domain&), GridView> indicatorFunction_;
+ // static const auto indicatorFunction_;
+
// VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
// VectorCT phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
// FieldVector<double,3> m_1_, m_2_, m_3_; // Corrector m_i coefficient vectors
@@ -105,8 +122,23 @@ public:
std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material");
Python::Module module = Python::import(materialFunctionName_);
- elasticityTensor_ = Python::make_function<MatrixRT>(module.get("H"));
-
+ elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
+
+ // module.get("Phases").toC<int>(Phases_);
+
+ auto indicatorFunction = Python::make_function<double>(module.get("indicatorFunction"));
+ indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
+
+
+
+ L1_ = Python::make_function<MatrixRT>(module.get("L1"));
+ L2_ = Python::make_function<MatrixRT>(module.get("L2"));
+ L3_ = Python::make_function<MatrixRT>(module.get("L3"));
+
+ // indicatorFunction_ = localFunction(indicatorFunctionGVF);
+
+
+
// Func2TensorParam elasticityTensor_ = Python::make_function<double>(module.get("L"));
// Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
// bool isotropic_ = true; // read from module File TODO
@@ -117,13 +149,27 @@ public:
- MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& 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);
+
+ }
+
// -----------------------------------------------------------------
@@ -148,6 +194,11 @@ public:
Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
+
+ // auto getIndicatorFunction() const {return indicatorFunction_;}
+ auto getIndicatorFunction() const {return localFunction(indicatorFunction_);}
+
+
// shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}
diff --git a/geometries/material.py b/geometries/material.py
index 4aba62b8..ff9bcfe2 100644
--- a/geometries/material.py
+++ b/geometries/material.py
@@ -1,16 +1,29 @@
import math
-
+from python_matrix_operations import *
+import ctypes
+import os
+import sys
Phases = 3
-lu = [1,2,3]
-mu_ = [3, 5, 6]
-lambda_ = [9, 7, 8]
+mu_ = [80, 80, 60]
+lambda_ = [80, 80, 25]
+
+# print('mu_:', mu_)
+# A = [[1, 5, 0], [5,1,0], [5,0,1]]
+#
+# print("sym(A):", sym(A))
-# lu = mu[0] mu[1] mu[2]
+# dir_path = os.path.dirname(os.path.realpath("/home/klaus/Desktop/Dune-Testing/dune-microstructure/dune/microstructure/matrix_operations.hh"))
+# handle = ctypes.CDLL(dir_path)
+#
+# handle.create2Darray.argtypes = [ctypes.c_int, ctypes.c_double, ctypes.c_double]
+#
+# def create2Darray(nside, mx, my):
+# return handle.create2Darray(nside, mx, my)
#Indicator function that determines both phases
@@ -21,23 +34,54 @@ lambda_ = [9, 7, 8]
###############
# Wood
###############
-def f(x):
+# def f(x):
+# theta=0.25
+# # mu_ = [3, 5, 6]
+# # lambda_ = [9, 7, 8]
+# # mu_ = 3 5 6
+# # lambda_ = 9 7 8
+
+# if ((abs(x[0]) < theta/2) and x[2]<0.25):
+# return [mu_[0], lambda_[0]] #latewood
+# # return 5 #latewood
+# elif ((abs(x[0]) > theta/2) and x[2]<0.25):
+# return [mu_[1], lambda_[1]] #latewood
+# # return 2
+# else :
+# return [mu_[2],lambda_[2]] #latewood #Phase3
+# # return 1
+
+def indicatorFunction(x):
theta=0.25
- # mu_ = [3, 5, 6]
- # lambda_ = [9, 7, 8]
- # mu_ = 3 5 6
- # lambda_ = 9 7 8
-
- if ((abs(x[0]) < theta/2) and x[2]<0.25):
- return [mu_[0], lambda_[0]] #latewood
- # return 5 #latewood
- elif ((abs(x[0]) > theta/2) and x[2]<0.25):
- return [mu_[1], lambda_[1]] #latewood
- # return 2
+ factor=1
+ if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+ return 1 #Phase1
+ elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+ return 2 #Phase2
else :
- return [mu_[2],lambda_[2]] #latewood #Phase3
- # return 1
+ return 0 #Phase3
+
+
+def L1(G):
+ return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
+def L2(G):
+ return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
+
+def L3(G):
+ return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
+
+
+# TEST
+
+# def L1(G):
+# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) )) #Phase1
+
+# def L2(G):
+# return Add(smult(2.0 * mu_[1], sym(G)),smult(lambda_[1] ,smult(trace(sym(G)),Id()) )) #Phase1
+
+# def L3(G):
+# return Add(smult(2.0 * mu_[2], sym(G)),smult(lambda_[2] ,smult(trace(sym(G)),Id()) )) #Phase1
#Workaround
@@ -47,14 +91,6 @@ def muValue(x):
def lambdaValue(x):
return lambda_
-# def b1(x):
-# return [[.5, 0, 0], [0,1,0], [0,0,0]]
-
-# def b2(x):
-# return [[.4, 0, 0], [0,.4,0], [0,0,0]]
-
-# def b3(x):
-# return [[0, 0, 0], [0,0,0], [0,0,0]]
@@ -99,15 +135,99 @@ def b3(x):
# --- elasticity tensor
# def L(G,x):
+# def L(G):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+
+
+
+
+
+
+
-def L(G):
+# --- elasticity tensor
+def L(G,x):
# input: symmetric matrix G, position x
# output: symmetric matrix LG
- return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+ theta=0.25
+ if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+ return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
+ elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+ return 2.0 * mu_[1] * sym(G) + lambda_[1] * trace(sym(G)) * Id() #Phase2
+ else :
+ return 2.0 * mu_[2] * sym(G) + lambda_[2] * trace(sym(G)) * Id() #Phase3
+# # # 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
+
+
+# def L(G,x):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# theta=0.25
+# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+# # return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
+# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) ))
+# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) )) #Phase2
+# else :
+# return Add(smult(2.0 * mu_[0], sym(G)),smult(lambda_[0] ,smult(trace(sym(G)),Id()) )) #Phase3
+# # return [[0, 0, 0], [0,0,0], [0,0,0]] #Phase3
+
+
+##TEST
+# def L(G,x):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# theta=0.25
+# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+# return [[x[0], 1, x[0]], [1, 1, 1],[x[0], x[0], x[0]]]
+# else :
+# return [[0, x[2], x[2]], [0,x[2],0], [0,0,0]]
+
+
+##TEST
+# def L(G,x):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# theta=0.25
+# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+# return sym([[1, 1, 1], [1, 1, 1],[1, 1, 1]])
+# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+# return sym([[x[0], 1, x[0]], [1, 1, 1],[x[0], x[0], x[0]]])
+# else :
+# return sym([[0, x[2], x[2]], [0,x[2],0], [0,0,0]])
+# # small speedup..
+# def L(G,x):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# theta=0.25
+# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+# return mu_[0] * (np.array(G).transpose() + np.array(G)) + lambda_[0] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3) #Phase1
+# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+# return mu_[1] * (np.array(G).transpose() + np.array(G)) + lambda_[1] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3) #Phase2
+# else :
+# return mu_[2] * (np.array(G).transpose() + np.array(G)) + lambda_[2] * (G[0][0] + G[1][1] + G[2][2]) * np.identity(3) #Phase3
+# # 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
+
+
+
+
+
+
+# def H(G,x):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# if (abs(x[0]) > 0.25):
+# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+# else:
+# return [[0, 0, 0], [0,0,0], [0,0,0]]
def H(G,x):
@@ -117,3 +237,5 @@ def H(G,x):
return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
else:
return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
diff --git a/geometries/python_matrix_operations.py b/geometries/python_matrix_operations.py
new file mode 100644
index 00000000..e160f3b3
--- /dev/null
+++ b/geometries/python_matrix_operations.py
@@ -0,0 +1,46 @@
+import numpy as np
+
+
+def sym(A): # 1/2 (A^T + A)
+ return 0.5 * (np.array(A).transpose() + np.array(A) )
+
+
+def trace(A):
+ return A[0][0] + A[1][1] + A[2][2]
+
+def Id():
+ return np.identity(3)
+
+def mult(A,B):
+ tmp = [[0, 0, 0], [0,0,0], [0,0,0]]
+ # iterate through rows of X
+ for i in range(3):
+ # iterate through columns of Y
+ for j in range(3):
+ # iterate through rows of Y
+ for k in range(3):
+ tmp[i][j] += A[i][k] * B[k][j]
+ return tmp
+
+def Add(A,B):
+ tmp = [[0, 0, 0], [0,0,0], [0,0,0]]
+ for i in range(3):
+ for j in range(3):
+ tmp[i][j] = (A[i][j] + B[i][j])
+ return tmp
+
+
+def smult(k,A):
+ return [[k*A[0][0], k*A[0][1], k*A[0][2]], [k*A[1][0],k*A[1][1],k*A[1][2]], [k*A[2][0],k*A[2][1],k*A[2][2]]]
+
+# def Id():
+# return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+
+
+# def sym(A): # 1/2 (A^T + A)
+# tmp = [[0, 0, 0], [0,0,0], [0,0,0]]
+# for i in range(3):
+# for j in range(3):
+# tmp[i][j] = 0.5 *(A[i][j] + A[j][i])
+# return tmp
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 6e4b2cd6..72a3744f 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -58,6 +58,8 @@
#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>
@@ -111,8 +113,8 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
-// a function:
-int half(int x, int y) {return x/2+y/2;}
+// // a function:
+// int half(int x, int y) {return x/2+y/2;}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -310,7 +312,7 @@ int main(int argc, char *argv[])
// 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"));
+ // auto materialFunction_ = Python::make_function<FieldVector<double,2>>(module.get("f"));
int Phases;
module.get("Phases").toC<int>(Phases);
@@ -341,19 +343,35 @@ int main(int argc, char *argv[])
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&) >;
+ 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"));
- // MatrixDomainFunc TestTensor = Python::make_function<MatrixRT>(module.get("H"));
// auto elasticityTensorGVF = Dune::Functions::makeGridViewFunction(elasticityTensor , Basis_CE.gridView());
// auto localElasticityTensor = localFunction(elasticityTensorGVF);
@@ -366,6 +384,9 @@ int main(int argc, char *argv[])
// 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 );
@@ -378,19 +399,23 @@ int main(int argc, char *argv[])
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));
- // printmatrix(std::cout, temp2, "temp2", "--");
-
+ // printmatrix(std::cout, tmp3, "tmp3", "--");
}
}
@@ -407,8 +432,8 @@ int main(int argc, char *argv[])
- std::function<int(int,int)> fn1 = half;
- std::cout << "fn1(60,20): " << fn1(60,20) << '\n';
+ // std::function<int(int,int)> fn1 = half;
+ // std::cout << "fn1(60,20): " << fn1(60,20) << '\n';
// std::cout << typeid(elasticityTensorGVF).name() << '\n';
@@ -444,27 +469,27 @@ int main(int argc, char *argv[])
//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 << "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)
@@ -531,12 +556,17 @@ int main(int argc, char *argv[])
// };
-
//------------------------------------------------------------------------------------------------
//--- compute Correctors
- auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
+ // 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();
@@ -549,10 +579,14 @@ int main(int argc, char *argv[])
if(print_debug)
correctorComputer.checkSymmetry();
+
//--- compute effective quantities
- auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,B_Term);
+ 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
@@ -748,4 +782,7 @@ int main(int argc, char *argv[])
log.close();
std::cout << "Total time elapsed: " << globalTimer.elapsed() << std::endl;
+
+
+
}
diff --git a/src/deprecated_code/elasticityTensor-globalFunctionVersion/Cell-Problem-New.cc b/src/deprecated_code/elasticityTensor-globalFunctionVersion/Cell-Problem-New.cc
new file mode 100644
index 00000000..09859eb5
--- /dev/null
+++ b/src/deprecated_code/elasticityTensor-globalFunctionVersion/Cell-Problem-New.cc
@@ -0,0 +1,766 @@
+#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 <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>(); //NEW
+ 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"));
+
+ // 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()))
+ {
+
+ 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 << "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));
+ // 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();
+
+
+ //--- 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/deprecated_code/elasticityTensor-globalFunctionVersion/CorrectorComputer.hh b/src/deprecated_code/elasticityTensor-globalFunctionVersion/CorrectorComputer.hh
new file mode 100644
index 00000000..633b2e7f
--- /dev/null
+++ b/src/deprecated_code/elasticityTensor-globalFunctionVersion/CorrectorComputer.hh
@@ -0,0 +1,1325 @@
+#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 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_);}
+
+
+ // --- 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();
+
+ // 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));
+ // printmatrix(std::cout, etmp, "etmp", "--");
+ double energyDensity= scalarProduct(etmp,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++)
+ {
+ // double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
+ double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(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;
+
+ 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 = 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>;
+
+ // 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_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(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_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),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 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);
+
+ 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 = 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! apparently (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/elasticityTensor-globalFunctionVersion/EffectiveQuantitiesComputer.hh b/src/deprecated_code/elasticityTensor-globalFunctionVersion/EffectiveQuantitiesComputer.hh
new file mode 100644
index 00000000..6c9d38ce
--- /dev/null
+++ b/src/deprecated_code/elasticityTensor-globalFunctionVersion/EffectiveQuantitiesComputer.hh
@@ -0,0 +1,469 @@
+#ifndef DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
+#define DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
+
+#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>
+#include <dune/common/parametertree.hh>
+
+using namespace Dune;
+using namespace MatrixOperations;
+using std::shared_ptr;
+using std::make_shared;
+using std::string;
+using std::cout;
+using std::endl;
+
+// template <class Basis>
+// class EffectiveQuantitiesComputer : public CorrectorComputer<Basis,Material> {
+
+template <class Basis, class Material>
+class EffectiveQuantitiesComputer {
+
+public:
+ static const int dimworld = 3;
+ // static const int nCells = 4;
+
+ 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;
+
+ using Func2Tensor = typename CorrectorComputer<Basis,Material>::Func2Tensor;
+ using FuncVector = typename CorrectorComputer<Basis,Material>::FuncVector;
+
+ using VectorCT = typename CorrectorComputer<Basis,Material>::VectorCT;
+
+ using HierarchicVectorView = typename CorrectorComputer<Basis,Material>::HierarchicVectorView;
+
+protected:
+
+ CorrectorComputer<Basis,Material>& correctorComputer_;
+ Func2Tensor prestrain_;
+ const Material& material_;
+
+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
+
+
+ // 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_;
+
+ // 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)
+ {
+
+ // 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
+ ///////////////////////////////
+ CorrectorComputer<Basis,Material> getCorrectorComputer(){return correctorComputer_;}
+
+ const shared_ptr<Basis> getBasis()
+ {
+ return correctorComputer_.getBasis();
+ }
+
+ auto getQeff(){return Q_;}
+ auto getBeff(){return Beff_;}
+
+
+ // -----------------------------------------------------------------
+ // --- Compute Effective Quantities
+ void computeEffectiveQuantities()
+ {
+
+ // Get everything.. better TODO: with Inheritance?
+ // auto test = correctorComputer_.getLoad_alpha1();
+ // auto phiContainer = correctorComputer_.getPhicontainer();
+ auto MContainer = correctorComputer_.getMcontainer();
+ auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
+ auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
+ auto mu_ = *correctorComputer_.getMu();
+ auto lambda_ = *correctorComputer_.getLambda();
+ auto gamma = correctorComputer_.getGamma();
+ auto basis = *correctorComputer_.getBasis();
+ ParameterTree parameterSet = correctorComputer_.getParameterSet();
+
+ shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(),
+ correctorComputer_.getCorr_phi2(),
+ correctorComputer_.getCorr_phi3()
+ };
+
+ auto prestrainGVF = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
+ auto prestrainFunctional = localFunction(prestrainGVF);
+
+ Q_ = 0 ;
+ Bhat_ = 0;
+
+ for(size_t a=0; a < 3; a++)
+ for(size_t b=0; b < 3; b++)
+ {
+ double energy = 0.0;
+ double prestrain = 0.0;
+ auto localView = basis.localView();
+ // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
+ auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
+ // auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,phiContainer[b]));
+ auto localfun_a = localFunction(GVFunc_a);
+ // auto localfun_b = localFunction(GVFunc_b);
+
+ ///////////////////////////////////////////////////////////////////////////////
+
+ auto matrixFieldG1GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[a], basis.gridView());
+ auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
+ 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);
+
+ // using GridView = typename Basis::GridView;
+
+ for (const auto& e : elements(basis.gridView()))
+ {
+ localView.bind(e);
+ matrixFieldG1.bind(e);
+ matrixFieldG2.bind(e);
+ localfun_a.bind(e);
+ // DerPhi2.bind(e);
+ mu.bind(e);
+ lambda.bind(e);
+ prestrainFunctional.bind(e);
+
+ double elementEnergy = 0.0;
+ double elementPrestrain = 0.0;
+
+ 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);
+ // 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 Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a];
+
+
+ auto G1 = matrixFieldG1(quadPos);
+ auto G2 = matrixFieldG2(quadPos);
+ // auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
+ // auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
+
+ auto X1 = G1 + Chi1;
+ // auto X2 = G2 + Chi2;
+
+
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
+ if (b==0)
+ {
+ elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos)) * quadPoint.weight() * integrationElement;
+ }
+ }
+ energy += elementEnergy;
+ prestrain += elementPrestrain;
+
+ }
+ Q_[a][b] = energy;
+ if (b==0)
+ Bhat_[a] = prestrain;
+ }
+ if(parameterSet.get<bool>("print_debug", false))
+ {
+ printmatrix(std::cout, Q_, "Matrix Q_", "--");
+ printvector(std::cout, Bhat_, "Bhat_", "--");
+ }
+
+ ///////////////////////////////
+ // Compute effective Prestrain B_eff (by solving linear system)
+ //////////////////////////////
+
+ // std::cout << "------- Information about Q matrix -----" << std::endl; // TODO
+ // MatrixInfo<MatrixRT> matrixInfo(Q_,true,2,1);
+ // std::cout << "----------------------------------------" << std::endl;
+ Q_.solve(Beff_,Bhat_);
+ if(parameterSet.get<bool>("print_debug", false))
+ printvector(std::cout, Beff_, "Beff_", "--");
+
+
+ //LOG-Output
+ auto& log = *(correctorComputer_.getLog());
+ log << "--- Prestrain Output --- " << std::endl;
+ log << "Bhat_: " << Bhat_ << std::endl;
+ log << "Beff_: " << Beff_ << " (Effective Prestrain)" << std::endl;
+ log << "------------------------ " << std::endl;
+
+ // TEST
+ // std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
+ return ;
+ }
+
+
+ // -----------------------------------------------------------------
+ // --- write Data to Matlab / Optimization-Code
+ void writeToMatlab(std::string outputPath)
+ {
+ std::cout << "write effective quantities to Matlab folder..." << std::endl;
+ //writeMatrixToMatlab(Q, "../../Matlab-Programs/QMatrix.txt");
+ writeMatrixToMatlab(Q_, outputPath + "/QMatrix.txt");
+ // write effective Prestrain in Matrix for Output
+ FieldMatrix<double,1,3> BeffMat;
+ BeffMat[0] = Beff_;
+ writeMatrixToMatlab(BeffMat, outputPath + "/BMatrix.txt");
+ return;
+ }
+
+
+
+
+ template<class MatrixFunction>
+ double energySP(const MatrixFunction& matrixFieldFuncA,
+ const MatrixFunction& matrixFieldFuncB)
+ {
+ double energy = 0.0;
+ auto mu_ = *correctorComputer_.getMu();
+ auto lambda_ = *correctorComputer_.getLambda();
+ auto gamma = correctorComputer_.getGamma();
+ auto basis = *correctorComputer_.getBasis();
+ auto localView = basis.localView();
+
+ 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);
+ for (const auto& e : elements(basis.gridView()))
+ {
+ localView.bind(e);
+ matrixFieldA.bind(e);
+ matrixFieldB.bind(e);
+ mu.bind(e);
+ lambda.bind(e);
+
+ double elementEnergy = 0.0;
+
+ auto geometry = e.geometry();
+ 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(e.type(), orderQR);
+ for (const auto& quadPoint : quad)
+ {
+ const auto& quadPos = quadPoint.position();
+ const double integrationElement = geometry.integrationElement(quadPos);
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), matrixFieldA(quadPos), matrixFieldB(quadPos));
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+ }
+ energy += elementEnergy;
+ }
+ return energy;
+ }
+
+
+
+ // --- Alternative that does not use orthogonality relation (75) in the paper
+ // void computeFullQ()
+ // {
+ // auto MContainer = correctorComputer_.getMcontainer();
+ // auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
+ // auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
+ // auto mu_ = *correctorComputer_.getMu();
+ // auto lambda_ = *correctorComputer_.getLambda();
+ // auto gamma = correctorComputer_.getGamma();
+ // auto basis = *correctorComputer_.getBasis();
+
+ // shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(),
+ // correctorComputer_.getCorr_phi2(),
+ // correctorComputer_.getCorr_phi3()
+ // };
+
+ // auto prestrainGVF = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
+ // auto prestrainFunctional = localFunction(prestrainGVF);
+
+ // Q_ = 0 ;
+ // Bhat_ = 0;
+
+ // for(size_t a=0; a < 3; a++)
+ // for(size_t b=0; b < 3; b++)
+ // {
+ // double energy = 0.0;
+ // double prestrain = 0.0;
+ // auto localView = basis.localView();
+ // // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
+ // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
+ // auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[b]));
+ // auto localfun_a = localFunction(GVFunc_a);
+ // auto localfun_b = localFunction(GVFunc_b);
+
+ // ///////////////////////////////////////////////////////////////////////////////
+ // auto matrixFieldG1GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[a], basis.gridView());
+ // auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
+ // 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);
+
+ // // using GridView = typename Basis::GridView;
+
+ // for (const auto& e : elements(basis.gridView()))
+ // {
+ // localView.bind(e);
+ // matrixFieldG1.bind(e);
+ // matrixFieldG2.bind(e);
+ // localfun_a.bind(e);
+ // localfun_b.bind(e);
+ // mu.bind(e);
+ // lambda.bind(e);
+ // prestrainFunctional.bind(e);
+
+ // double elementEnergy = 0.0;
+ // double elementPrestrain = 0.0;
+
+ // 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);
+ // // 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 Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a] + matrixFieldG1(quadPos);
+ // auto Chi2 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_b(quadPos))) + *MContainer[b] + matrixFieldG2(quadPos);
+
+ // // auto G1 = matrixFieldG1(quadPos);
+ // // auto G2 = matrixFieldG2(quadPos);
+ // // auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
+ // // auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
+
+ // // auto X1 = G1 + Chi1;
+ // // auto X2 = G2 + Chi2;
+
+
+ // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), Chi1, Chi2);
+ // elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
+ // }
+ // energy += elementEnergy;
+ // prestrain += elementPrestrain;
+
+ // }
+ // Q_[a][b] = energy;
+ // if (b==0)
+ // Bhat_[a] = prestrain;
+ // }
+ // printmatrix(std::cout, Q_, "Matrix Q_", "--");
+ // printvector(std::cout, Bhat_, "Bhat_", "--");
+ // ///////////////////////////////
+ // // Compute effective Prestrain B_eff (by solving linear system)
+ // //////////////////////////////
+ // // std::cout << "------- Information about Q matrix -----" << std::endl; // TODO
+ // // MatrixInfo<MatrixRT> matrixInfo(Q_,true,2,1);
+ // // std::cout << "----------------------------------------" << std::endl;
+ // Q_.solve(Beff_,Bhat_);
+ // printvector(std::cout, Beff_, "Beff_", "--");
+
+ // //LOG-Output
+ // auto& log = *(correctorComputer_.getLog());
+ // log << "--- Prestrain Output --- " << std::endl;
+ // log << "Bhat_: " << Bhat_ << std::endl;
+ // log << "Beff_: " << Beff_ << " (Effective Prestrain)" << std::endl;
+ // log << "------------------------ " << std::endl;
+ // return ;
+ // }
+
+
+}; // end class
+
+
+
+#endif
\ No newline at end of file
diff --git a/src/deprecated_code/elasticityTensor-globalFunctionVersion/material.py b/src/deprecated_code/elasticityTensor-globalFunctionVersion/material.py
new file mode 100644
index 00000000..2487cfce
--- /dev/null
+++ b/src/deprecated_code/elasticityTensor-globalFunctionVersion/material.py
@@ -0,0 +1,138 @@
+import math
+from python_matrix_operations import *
+
+
+
+
+Phases = 3
+
+mu_ = [3, 5, 0]
+lambda_ = [9, 7, 0]
+
+
+
+#Indicator function that determines both phases
+# x[0] : y1-component
+# x[1] : y2-component
+# x[2] : x3-component
+# --- replace with your definition of indicatorFunction:
+###############
+# Wood
+###############
+def f(x):
+ theta=0.25
+ # mu_ = [3, 5, 6]
+ # lambda_ = [9, 7, 8]
+ # mu_ = 3 5 6
+ # lambda_ = 9 7 8
+
+ if ((abs(x[0]) < theta/2) and x[2]<0.25):
+ return [mu_[0], lambda_[0]] #latewood
+ # return 5 #latewood
+ elif ((abs(x[0]) > theta/2) and x[2]<0.25):
+ return [mu_[1], lambda_[1]] #latewood
+ # return 2
+ else :
+ return [mu_[2],lambda_[2]] #latewood #Phase3
+ # return 1
+
+
+
+#Workaround
+def muValue(x):
+ return mu_
+
+def lambdaValue(x):
+ return lambda_
+
+
+
+
+
+###############
+# Cross
+###############
+# def f(x):
+# theta=0.25
+# factor=1
+# if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+# return 1 #Phase1
+# elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+# return 2 #Phase2
+# else :
+# return 0 #Phase3
+
+
+
+
+
+# def f(x):
+# # --- replace with your definition of indicatorFunction:
+# if (abs(x[0]) > 0.25):
+# return 1 #Phase1
+# else :
+# return 0 #Phase2
+
+def b1(x):
+ return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def b2(x):
+ return [[1, 0, 0], [0,1,0], [0,0,1]]
+
+def b3(x):
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# mu=80 80 60
+
+# lambda=80 80 25
+
+
+# --- elasticity tensor
+# def L(G,x):
+# def L(G):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+
+
+
+
+
+
+
+
+# --- elasticity tensor
+def L(G,x):
+ # input: symmetric matrix G, position x
+ # output: symmetric matrix LG
+ theta=0.25
+ if (x[0] <-1/2+theta and x[2]<-1/2+theta):
+ return 2.0 * mu_[0] * sym(G) + lambda_[0] * trace(sym(G)) * Id() #Phase1
+ elif (x[1]< -1/2+theta and x[2]>1/2-theta):
+ return 2.0 * mu_[1] * sym(G) + lambda_[1] * trace(sym(G)) * Id() #Phase2
+ else :
+ return 2.0 * mu_[2] * sym(G) + lambda_[2] * trace(sym(G)) * Id() #Phase3
+ # 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
+
+
+
+
+
+# def H(G,x):
+# # input: symmetric matrix G, position x
+# # output: symmetric matrix LG
+# if (abs(x[0]) > 0.25):
+# return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+# else:
+# return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+
+def H(G,x):
+ # input: symmetric matrix G, position x
+ # output: symmetric matrix LG
+ if (abs(x[0]) > 0.25):
+ return [[1, 1, 1], [1, 1, 1],[1, 1, 1]]
+ else:
+ return [[0, 0, 0], [0,0,0], [0,0,0]]
+
+# 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
\ No newline at end of file
diff --git a/src/deprecated_code/elasticityTensor-globalFunctionVersion/prestrainedMaterial.hh b/src/deprecated_code/elasticityTensor-globalFunctionVersion/prestrainedMaterial.hh
new file mode 100644
index 00000000..b0c68e3b
--- /dev/null
+++ b/src/deprecated_code/elasticityTensor-globalFunctionVersion/prestrainedMaterial.hh
@@ -0,0 +1,199 @@
+#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
+#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
+
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/yaspgrid.hh>
+#include <dune/microstructure/matrix_operations.hh>
+
+#include <dune/fufem/dunepython.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 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 CellGridType = YaspGrid< dim, EquidistantOffsetCoordinates< double, dim>>;
+ // using Domain = typename CellGridType::LeafGridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using ScalarRT = FieldVector< double, 1>;
+ using VectorRT = FieldVector< double, dimworld>;
+ using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
+ using FuncScalar = std::function< double(const Domain&) >;
+ using Func2Tensor = std::function< MatrixRT(const Domain&) >;
+ using Func2TensorParam = std::function< MatrixRT(const MatrixRT& ,const Domain&) >;
+
+
+protected:
+
+ const GridView& gridView_;
+ const ParameterTree& parameterSet_;
+
+
+ // const FieldVector<double , ...number of mu-Values/Phases> .. schwierig zur compile-time
+
+ // const FuncScalar mu_;
+ // const FuncScalar lambda_;
+ // double gamma_;
+
+ std::string materialFunctionName_;
+
+ // --- Number of material phases?
+ // const int phases_;
+
+ // Func2Tensor materialFunction_; //actually not needed??
+
+ // Func2Tensor elasticityTensor_;
+ Func2TensorParam elasticityTensor_;
+
+// 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_};
+
+
+public:
+ ///////////////////////////////
+ // constructor
+ ///////////////////////////////
+ prestrainedMaterial(const GridView gridView,
+ const ParameterTree& parameterSet) // string: "name of material"? // mu_(mu), muValues? müsste Anzahl Phasen bereits kennen..
+ : gridView_(gridView),
+ parameterSet_(parameterSet)
+ {
+ std::string materialFunctionName_ = parameterSet.get<std::string>("materialFunction", "material");
+ Python::Module module = Python::import(materialFunctionName_);
+
+ elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
+
+ // Func2TensorParam elasticityTensor_ = Python::make_function<double>(module.get("L"));
+ // Func2Tensor materialFunction_ = Python::make_function<double>(module.get("f"));
+ // bool isotropic_ = true; // read from module File TODO
+ // Func2Tensor elasticityTensor_ = Python::make_function<double>(module.get("L"));
+ // Func2Tensor elasticityTensor_ = Python::make_function<MatrixRT>(module.get("L"));
+ }
+
+
+
+
+ MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const
+ {
+ //--- apply elasticityTensor_ to input Matrix G at position x
+ return elasticityTensor_(G,x);
+
+ }
+
+
+
+ // -----------------------------------------------------------------
+ // --- 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_;}
+
+
+ // shared_ptr<Func2TensorParam> getElasticityTensor(){return make_shared<Func2TensorParam>(elasticityTensor_);}
+
+
+ //TODO getLameParameters() .. Throw Exception if isotropic_ = false!
+
+
+
+
+ // shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
+ // shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
+ // shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
+ // shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
+
+ // shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
+ // shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
+
+
+ // --- Get Correctors
+ // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
+ // auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
+ // auto getMcontainer(){return mContainer;}
+ // shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}
+
+
+ // // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
+ // auto getMatrixBasiscontainer(){return make_shared<std::array<MatrixRT,3 >>(MatrixBasisContainer_);}
+ // // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
+ // auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
+
+
+
+
+
+ // shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
+ // shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
+ // shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
+ // shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
+
+
+
+
+
+
+}; // end class
+
+
+
+
+#endif
--
GitLab