From 016d5c367d07f69c7ca141ef358b1e08b9ef871b Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Thu, 9 Sep 2021 18:52:21 +0200
Subject: [PATCH] Add Stripped-down version of Cell-Problem for the computation
of mu_gamma (significant speed-up)
---
src/CMakeLists.txt | 1 +
src/Cell-Problem_muGamma.cc | 1142 +++++++++++++++++++++++++++++++++++
2 files changed, 1143 insertions(+)
create mode 100644 src/Cell-Problem_muGamma.cc
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 402c96f0..73e5b2f8 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -4,6 +4,7 @@
set(programs Cell-Problem
+ Cell-Problem_muGamma
)
diff --git a/src/Cell-Problem_muGamma.cc b/src/Cell-Problem_muGamma.cc
new file mode 100644
index 00000000..5c620bdd
--- /dev/null
+++ b/src/Cell-Problem_muGamma.cc
@@ -0,0 +1,1142 @@
+#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/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/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/functions/functionspacebases/hierarchicvectorwrapper.hh>
+
+// #include <dune/fufem/discretizationerror.hh>
+// #include <boost/multiprecision/cpp_dec_float.hpp>
+// #include <boost/any.hpp>
+
+#include <any>
+#include <variant>
+#include <string>
+#include <iomanip>
+
+using namespace Dune;
+using namespace MatrixOperations;
+
+// ------------------------------------------------------------------------
+//
+// This is a stripped-down Version of Cell-Problem.cc
+// that only computes mu_gamma = q_3 in a faster manner
+// i.e. only the third corrector phi_3 is needed.
+// ------------------------------------------------------------------------
+
+
+
+
+
+
+
+
+
+
+
+
+
+//////////////////////////////////////////////////////////////////////
+// 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();
+}
+
+
+
+
+
+template<class Basis>
+auto arbitraryComponentwiseIndices(const Basis& basis,
+ 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;
+}
+
+
+
+
+template<class LocalView, class Matrix, class localFunction1, class localFunction2>
+void computeElementStiffnessMatrix(const LocalView& localView,
+ Matrix& elementMatrix,
+ const localFunction1& mu,
+ const localFunction2& lambda,
+ const double gamma
+ )
+{
+ // Local StiffnessMatrix of the form:
+ // | phi*phi m*phi |
+ // | phi *m m*m |
+ using Element = typename LocalView::Element;
+ const Element element = localView.element();
+ auto geometry = element.geometry();
+ constexpr int dim = Element::dimension;
+ constexpr int dimWorld = dim;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+
+ elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
+ elementMatrix = 0;
+
+ // 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}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+ for (const auto& quadPoint : quad)
+ {
+ 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 k=0; k < dimWorld; k++)
+ for (size_t i=0; i < nSf; i++)
+ {
+ // "phi*phi"-part
+ for (size_t l=0; l< dimWorld; l++)
+ for (size_t j=0; j < nSf; j++ )
+ {
+ // (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
+ // printmatrix(std::cout, defGradientU , "defGradientU", "--");
+ // (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
+
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
+// double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
+
+ size_t col = localView.tree().child(k).localIndex(i);
+ size_t row = localView.tree().child(l).localIndex(j);
+ elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
+
+ // "m*phi" & "phi*m" - part
+ for( size_t m=0; m<3; m++)
+ {
+ double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
+
+ 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++)
+ for(size_t n=0; n<3; n++)
+ {
+ double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
+ elementMatrix[localPhiOffset+m][localPhiOffset+n]= energyDensityGG * quadPoint.weight() * integrationElement;
+ }
+ }
+}
+
+
+
+// Get the occupation pattern of the stiffness matrix
+template<class Basis, class ParameterSet>
+void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet& parameterSet)
+{
+ // 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.template get<bool>("set_oneBasisFunction_Zero ", true)){
+ FieldVector<int,3> row;
+ unsigned int arbitraryLeafIndex = parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
+ unsigned int arbitraryElementNumber = parameterSet.template 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(basis,arbitraryElementNumber,arbitraryLeafIndex);
+
+ for (int k = 0; k<3; k++)
+ nb.add(row[k],row[k]);
+ }
+}
+
+
+
+// Compute the source term for a single element
+// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
+template<class LocalView, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void computeElementLoadVector( const LocalView& localView,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
+ const double gamma,
+ Vector& elementRhs,
+ const Force& 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}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ 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); // für simplex
+ const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+ 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]);
+
+ // "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]; // X3
+
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientV, forceTerm(quadPos));
+ 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 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], forceTerm(quadPos));
+ elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
+ }
+ }
+}
+
+
+
+template<class Basis, class Matrix, class LocalFunction1, class LocalFunction2, class ParameterSet>
+void assembleCellStiffness(const Basis& basis,
+ LocalFunction1& muLocal,
+ LocalFunction2& lambdaLocal,
+ const double gamma,
+ Matrix& matrix,
+ ParameterSet& parameterSet
+ )
+{
+ std::cout << "assemble Stiffness-Matrix begins." << std::endl;
+
+ MatrixIndexSet occupationPattern;
+ getOccupationPattern(basis, occupationPattern, parameterSet);
+ occupationPattern.exportIdx(matrix);
+ matrix = 0;
+
+ auto localView = basis.localView();
+ const int phiOffset = basis.dimension();
+
+ for (const auto& element : elements(basis.gridView()))
+ {
+ localView.bind(element);
+ muLocal.bind(element);
+ lambdaLocal.bind(element);
+ const int localPhiOffset = localView.size();
+
+ //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+ Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
+ computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal, gamma);
+ //printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
+ //std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
+ //std::cout << "elementMatrix.M() : " << elementMatrix.M() << 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", "--");
+ }
+}
+
+
+template<class Basis, class LocalFunction1, class LocalFunction2, class Vector, class Force>
+void assembleCellLoad(const Basis& basis,
+ LocalFunction1& muLocal,
+ LocalFunction2& lambdaLocal,
+ const double gamma,
+ Vector& b,
+ const Force& 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);
+
+ 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;
+
+ BlockVector<FieldVector<double,1> > elementRhs;
+ computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, loadFunctional);
+// 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", "--");
+ }
+}
+
+
+
+template<class Basis, class LocalFunction1, class LocalFunction2, class MatrixFunction>
+auto energy(const Basis& basis,
+ LocalFunction1& mu,
+ LocalFunction2& lambda,
+ const MatrixFunction& matrixFieldFuncA,
+ const MatrixFunction& matrixFieldFuncB)
+{
+
+// TEST HIGHER PRECISION
+// using float_50 = boost::multiprecision::cpp_dec_float_50;
+// float_50 higherPrecEnergy = 0.0;
+ double energy = 0.0;
+ constexpr int dim = Basis::LocalView::Element::dimension;
+ constexpr int dimWorld = dim;
+
+ 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);
+
+ using GridView = typename Basis::GridView;
+ using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+// 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);
+ matrixFieldA.bind(e);
+ matrixFieldB.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 + 5 ); // TEST
+ 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);
+
+ auto strain1 = matrixFieldA(quadPos);
+ auto strain2 = matrixFieldB(quadPos);
+
+ double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), strain1, strain2);
+
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+ //elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
+ }
+ energy += elementEnergy;
+ //higherPrecEnergy += elementEnergy_HP;
+ }
+// TEST
+// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
+ return energy;
+}
+
+
+
+template<class Basis, class Matrix, class Vector, class ParameterSet> // changed to take only one load...
+void setOneBaseFunctionToZero(const Basis& basis,
+ Matrix& stiffnessMatrix,
+ Vector& load_alpha,
+ ParameterSet& parameterSet
+ )
+{
+ 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(basis,arbitraryElementNumber,arbitraryLeafIndex);
+
+ for (int k = 0; k < dim; k++)
+ {
+ load_alpha[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;
+ }
+}
+
+
+
+template<class Basis>
+auto childToIndexMap(const Basis& basis,
+ 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;
+}
+
+
+template<class Basis, class Vector>
+auto integralMean(const Basis& basis,
+ Vector& coeffVector
+ )
+{
+ // computes integral mean of given LocalFunction
+
+ constexpr int dim = Basis::LocalView::Element::dimension;
+
+ 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);
+ 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 ;
+}
+
+
+template<class Basis, class Vector>
+auto subtractIntegralMean(const Basis& basis,
+ Vector& coeffVector
+ )
+{
+ // Substract correct Integral mean from each associated component function
+ auto IM = integralMean(basis, coeffVector);
+
+ constexpr int dim = Basis::LocalView::Element::dimension;
+
+ for(size_t k=0; k<dim; k++)
+ {
+ //std::cout << "Integral-Mean: " << IM[k] << std::endl;
+ auto idx = childToIndexMap(basis,k);
+
+ for ( int i : idx)
+ coeffVector[i] -= IM[k];
+ }
+}
+
+
+
+//////////////////////////////////////////////////
+// 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]))
+ );
+ };
+
+
+
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[])
+{
+ MPIHelper::instance(argc, argv);
+
+ 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/output.txt");
+// std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt");
+ std::string outputPath = parameterSet.get("outputPath", "/home/klaus/Desktop/DUNE/dune-microstructure/outputs");
+// std::string MatlabPath = parameterSet.get("MatlabPath", "/home/klaus/Desktop/DUNE/dune-microstructure/Matlab-Programs");
+// std::string outputPath = "/home/klaus/Desktop/DUNE/dune-microstructure/outputs/output.txt";
+ std::fstream log;
+ log.open(outputPath + "/output.txt" ,std::ios::out);
+
+ std::cout << "outputPath:" << outputPath << std::endl;
+
+// parameterSet.report(log); // short Alternativ
+
+ constexpr int dim = 3;
+ constexpr int dimWorld = 3;
+
+ ///////////////////////////////////
+ // Get Parameters/Data
+ ///////////////////////////////////
+ double gamma = parameterSet.get<double>("gamma",3.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");
+ double beta = parameterSet.get<double>("beta",2.0);
+ double mu1 = parameterSet.get<double>("mu1",10.0);;
+ double mu2 = beta*mu1;
+ double lambda1 = parameterSet.get<double>("lambda1",0.0);;
+ double lambda2 = beta*lambda1;
+
+ // Plate geometry settings
+ double width = parameterSet.get<double>("width", 1.0); //geometry cell, cross section
+ // double len = parameterSet.get<double>("len", 10.0); //length
+ // double height = parameterSet.get<double>("h", 0.1); //rod thickness
+ // double eps = parameterSet.get<double>("eps", 0.1); //size of perioticity cell
+
+ ///////////////////////////////////
+ // Get Prestrain/Parameters
+ ///////////////////////////////////
+// unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", "analytical_Example"); //OLD
+// std::string prestraintype = parameterSet.get<std::string>("prestrainType", "analytical_Example");
+// double rho1 = parameterSet.get<double>("rho1", 1.0);
+// double rho2 = alpha*rho1;
+// auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
+// auto B_Term = prestrainImp.getPrestrain(prestraintype);
+
+
+
+ log << "----- Input Parameters -----: " << std::endl;
+// log << "prestrain imp: " << prestraintype << "\nrho1 = " << rho1 << "\nrho2 = " << rho2 << 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:
+ unsigned int cellDomain = parameterSet.get<unsigned int>("cellDomain", 1);
+ // (shifted) 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});
+ if (cellDomain == 2)
+ {
+ // Domain : [0,1)^2 x (-1/2, 1/2)
+ FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
+ FieldVector<double,dim> upper({1.0, 1.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:: #1 level #2 |q1_a-q1_c| #3 |q2_a-q2_c| #4 |q3_a-q3_c| #5 |b1_a-b1_c| #6 |b2_a-b2_c| #7 |b3_a-b3_c|
+ 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 << "Level: " << 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 Elements in each direction: " << nElements << std::endl;
+ log << "Number of 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();
+
+ using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
+ using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
+ 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();
+ 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 Solver ---
+ // 1 : CG-Solver
+ // 2 : GMRES
+ // 3 : QR
+ ///////////////////////////////////
+ unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
+
+ // Print Options
+ bool print_debug = parameterSet.get<bool>("print_debug", false);
+
+
+ bool write_materialFunctions = parameterSet.get<bool>("write_materialFunctions", false);
+ bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false);
+
+
+ bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
+ bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
+ bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
+
+ /////////////////////////////////////////////////////////
+ // Choose a finite element space for Cell Problem
+ /////////////////////////////////////////////////////////
+ using namespace Functions::BasisFactory;
+ Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
+
+ // 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)});
+ }
+
+ // First order periodic Lagrange-Basis
+ auto Basis_CE = makeBasis(
+ gridView_CE,
+ power<dim>(
+ Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices),
+ flatLexicographic()));
+
+ log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
+ const int phiOffset = Basis_CE.size();
+
+ /////////////////////////////////////////////////////////
+ // Data structures: Stiffness matrix and right hand side vectors
+ /////////////////////////////////////////////////////////
+ VectorCT load_alpha1, load_alpha2, load_alpha3;
+ MatrixCT stiffnessMatrix_CE;
+
+ bool set_IntegralZero = parameterSet.get<bool>("set_IntegralZero", true);
+ bool set_oneBasisFunction_Zero = false;
+ bool substract_integralMean = false;
+ if(set_IntegralZero)
+ {
+ set_oneBasisFunction_Zero = true;
+ substract_integralMean = true;
+ }
+
+ /////////////////////////////////////////////////////////
+ // Define "rhs"-functions
+ /////////////////////////////////////////////////////////
+ Func2Tensor x3G_3 = [] (const Domain& x) {
+ return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ };
+
+ Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {return -1.0*x3G_3(x);};
+
+
+
+ ///////////////////////////////////////////////
+ // 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}};
+ MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
+ // printmatrix(std::cout, basisContainer[0] , "G_1", "--");
+ // printmatrix(std::cout, basisContainer[1] , "G_2", "--");
+ // printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
+ // log << basisContainer[0] << std::endl;
+
+
+ //////////////////////////////////////////////////////////////////////
+ // Determine global indices of components for arbitrary (local) index
+ //////////////////////////////////////////////////////////////////////
+ int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts
+ int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
+
+ FieldVector<int,3> row;
+ row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
+ if(print_debug)
+ printvector(std::cout, row, "row" , "--");
+
+ /////////////////////////////////////////////////////////
+ // Assemble the system
+ /////////////////////////////////////////////////////////
+ assembleCellStiffness(Basis_CE, muLocal, lambdaLocal, gamma, stiffnessMatrix_CE, parameterSet);
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg); // only third corrector is needed
+ // printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
+ // printvector(std::cout, load_alpha1, "load_alpha1", "--");
+
+
+
+ // set one basis-function to zero
+ if(set_oneBasisFunction_Zero)
+ {
+ setOneBaseFunctionToZero(Basis_CE, stiffnessMatrix_CE, load_alpha3, parameterSet);
+ }
+
+ ////////////////////////////////////////////////////
+ // Compute solution
+ ////////////////////////////////////////////////////
+ VectorCT x_3 = load_alpha3;
+
+ if (Solvertype == 1) // CG - SOLVER
+ {
+ std::cout << "------------ CG - Solver ------------" << std::endl;
+ MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_CE);
+
+
+
+ // Sequential incomplete LU decomposition as the preconditioner
+ SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_CE,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
+ 2); // 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;
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ 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_CE);
+
+ // 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
+ 2); // Verbosity of the solver
+
+ // Object storing some statistics about the solving process
+ InverseOperatorResult statistics;
+
+
+ solver.apply(x_3, load_alpha3, statistics);
+ log << "Solver-type used: " <<" GMRES-Solver" << 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_CE);
+ sPQR.setVerbosity(1);
+ InverseOperatorResult statisticsQR;
+
+ sPQR.apply(x_3, load_alpha3, statisticsQR);
+ log << "Solver-type used: " <<" QR-Solver" << std::endl;
+ }
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Extract phi_alpha & M_alpha coefficients
+ ////////////////////////////////////////////////////////////////////////////////////
+ VectorCT phi_3;
+ phi_3.resize(Basis_CE.size());
+ phi_3 = 0;
+
+ for(size_t i=0; i<Basis_CE.size(); i++)
+ {
+
+ phi_3[i] = x_3[i];
+ }
+
+ FieldVector<double,3> m_1, m_2, m_3;
+
+ for(size_t i=0; i<3; 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);
+
+ for(size_t i=0; i<3; i++)
+ {
+ M_3 += m_3[i]*basisContainer[i];
+ }
+
+ std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
+// printmatrix(std::cout, M_3, "Corrector-Matrix M_3", "--");
+ log << "---------- OUTPUT ----------" << std::endl;
+ log << " --------------------" << std::endl;
+ log << "Corrector-Matrix M_3: \n" << M_3 << std::endl;
+ log << " --------------------" << std::endl;
+
+ ////////////////////////////////////////////////////////////////////////////
+ // Substract Integral-mean
+ ////////////////////////////////////////////////////////////////////////////
+ using SolutionRange = FieldVector<double,dim>;
+
+ if(substract_integralMean)
+ {
+ std::cout << " --- substracting integralMean --- " << std::endl;
+ subtractIntegralMean(Basis_CE, phi_3);
+ subtractIntegralMean(Basis_CE, x_3);
+ }
+
+ ////////////////////////////////////////////////////////////////////////////
+ // Make a discrete function from the FE basis and the coefficient vector
+ //////////////////////////////////////////////////////////////////////////// //
+// auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
+// Basis_CE,
+// phi_3);
+
+
+ /////////////////////////////////////////////////////////
+ // Compute effective quantities: Elastic law & Prestrain
+ /////////////////////////////////////////////////////////
+ VectorCT tmp1;
+ assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3G_3);
+ double GGterm = 0.0;
+ GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3G_3, x3G_3 );
+ auto mu_gamma = (x_3*tmp1) + GGterm;
+ log << "mu_gamma=" << mu_gamma << std::endl;
+
+ std::cout << "mu_gamma:" << mu_gamma << std::endl;
+ Storage_Quantities.push_back(mu_gamma);
+
+
+ levelCounter++;
+ } // Level-Loop End
+
+
+ log.close();
+}
--
GitLab