Skip to content
Snippets Groups Projects
Commit 2d5ffc6e authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Add Stripped-down version of Cell-Problem for the computation of mu_gamma (significant speed-up)

parent 7e83133c
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/CMakeLists.txt
+ 1
− 0
View file @ 2d5ffc6e
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
set(programs Cell-Problem set(programs Cell-Problem
Cell-Problem_muGamma
) )
......
src/Cell-Problem_muGamma.cc 0 → 100644
+ 1142
− 0
View file @ 2d5ffc6e
#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();
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment