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

Add computation of MuGamma via the 2D Poisson-Type Equation

parent 3f8e7cdf
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/microstructure/matrix_operations.hh
+ 15
− 8
View file @ 3a86d4be
...@@ -83,14 +83,21 @@ namespace MatrixOperations { ...@@ -83,14 +83,21 @@ namespace MatrixOperations {
return sum; return sum;
} }
static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){ // static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){ // ?? Check with Robert
E1= sym(E1); // E1= sym(E1);
E2 = sym(E2); // E2 = sym(E2);
double t1 = scalarProduct(E1,E2); // double t1 = scalarProduct(E1,E2);
t1 *= 2* mu; // t1 *= 2* mu;
double t2 = trace(E1)*trace(E2); // double t2 = trace(E1)*trace(E2);
t2 *= lambda; // t2 *= lambda;
return t1 + t2; // return t1 + t2;
// }
static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
{
auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
return scalarProduct(t1,E2);
} }
// --- Generalization: Define Quadratic QuadraticForm // --- Generalization: Define Quadratic QuadraticForm
......
src/CMakeLists.txt
+ 1
− 0
View file @ 3a86d4be
...@@ -5,6 +5,7 @@ ...@@ -5,6 +5,7 @@
set(programs Cell-Problem set(programs Cell-Problem
Cell-Problem_muGamma Cell-Problem_muGamma
Compute_MuGamma
) )
......
src/Compute_MuGamma.cc 0 → 100644
+ 1043
− 0
View file @ 3a86d4be
#include <config.h>
#include <vector>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/gmshreader.hh>
#include <dune/grid/yaspgrid.hh>
#include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/io.hh>
#include <dune/istl/matrix.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/matrixmarket.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh> // use for LagrangeBasis
#include <dune/common/float_cmp.hh>
#include <dune/common/math.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/functions/functionspacebases/periodicbasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/functionspacebases/boundarydofs.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/microstructure/prestrain_material_geometry.hh>
// #include <dune/microstructure/matrix_operations.hh>
// #include <math.h>
#include <cmath>
// #include <numbers>
using namespace Dune;
// ------------------------------------------------------------------------
//
// Solving the 2D "Poisson-Type" Equation ( (38) in the draft)
// in order to compute mu_Gamma in a faster manner
//
// ------------------------------------------------------------------------
template<class LocalView, class Matrix, class LocalFunction>
void assembleElementStiffnessMatrix(const LocalView& localView,
Matrix& elementMatrix,
const LocalFunction& mu,
const double gamma
)
{
using Element = typename LocalView::Element;
constexpr int dim = Element::dimension;
const Element element = localView.element();
auto geometry = element.geometry();
const auto& localFiniteElement = localView.tree().finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// std::cout << "nSf:" << nSf << std::endl;
elementMatrix.setSize(localView.size(),localView.size());
elementMatrix = 0;
// int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1) + 5; // doesnt change anything
int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
const auto& quadRule = QuadratureRules<double, dim>::rule(element.type(),orderQR);
// std::cout << "QuadRule-Order:" << orderQR << std::endl;
// double muValue = 0.0;
// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
for (const auto& quadPoint : quadRule)
{
const auto& quadPos = quadPoint.position();
// std::cout << " quadPOS: " << quadPos << std::endl;
// TEST
// double muValue = mu(quadPos);
// std::cout << "muValue:" << muValue << std::endl;
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
const double 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++)
jacobian.mv(referenceGradients[i][0], gradients[i]); //TODO ? passt..
// // print all gradients //TEST
// FieldVector<double,dim> tmp = {0.0 , 0.0};
// for (size_t i=0; i < nSf; i++)
// {
// printvector(std::cout, gradients[i], "gradients[i]", "--" );
//
//
// tmp[0] += gradients[i][0];
// tmp[1] += gradients[i][1];
//
// // tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
// // tmp[1] += coeffVector[globalIdx]*gradients[i][2];
// // printvector(std::cout, tmp, "tmp", "--" );
//
// }
// printvector(std::cout, tmp, "sum of basisfunction Gradients", "--" );
for (size_t p=0; p<elementMatrix.N(); p++)
{
for (size_t q=0; q<elementMatrix.M(); q++)
{
// auto localRow = localView.tree().localIndex(p); // VERTAUSCHT?!?!
// auto localCol = localView.tree().localIndex(q);
auto localRow = localView.tree().localIndex(q);
auto localCol = localView.tree().localIndex(p);
// auto value = mu(quadPos)*(2.0* gradients[p][0] * gradients[q][0])+ mu(quadPos)*((1.0/(std::pow(gamma,2)))*(gradients[p][1] * gradients[q][1]));
// auto value = (mu(quadPos)*gradients[p][0] * gradients[q][0])+ ((1.0/gamma)*(gradients[p][1] * gradients[q][1]));
// std::cout << "gradients[p][0]" << gradients[p][0] << std::endl;
// std::cout << "gradients[q][0]" << gradients[q][0] << std::endl;
// std::cout << "(1.0/std::pow(gamma,2))*gradients[p][1]" << (1.0/std::pow(gamma,2))*gradients[p][1]<< std::endl;
// auto value3 = mu(quadPos)*(1.0/gamma)*gradients[p][2]; // 3D-Version
// auto value4 = value3*gradients[q][2]; // 3D-Version
auto value1 = 2.0*mu(quadPos)* gradients[p][0] *gradients[q][0]; //#1
// auto value1 = 2.0*mu(quadPos)*sqrt(2.0)* gradients[p][0] *gradients[q][0]; // TEST TODO warum passt es damit besser bei gamma = 1.0 ???
// auto value2 = mu(quadPos)*(1.0/std::pow(gamma,2))*gradients[p][1] * gradients[q][1] ;
auto value2 = 2.0*mu(quadPos)*(1.0/std::pow(gamma,2))*gradients[p][1] * gradients[q][1] ; //#1 TEST FAKTOR 2 hat hier gefehlt?!?!
auto value = value1 + value2;
elementMatrix[localRow][localCol] += value * quadPoint.weight() * integrationElement;
}
}
}
}
// Compute the source term for a single element
template<class LocalView,class Vector, class LocalFunction, class Force>
void assembleElementVolumeTerm(const LocalView& localView,
Vector& elementRhs,
LocalFunction& mu,
const Force& forceTerm)
{
using Element = typename LocalView::Element;
auto element = localView.element();
auto geometry = element.geometry();
constexpr int dim = Element::dimension;
// Set of shape functions for a single element
const auto& localFiniteElement = localView.tree().finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// Set all entries to zero
elementRhs.resize(localFiniteElement.size());
elementRhs = 0;
// A quadrature rule
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
// std::cout << "elementRhs.size():" << elementRhs.size() << std::endl;
const auto& quadRule = QuadratureRules<double,dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quadRule)
{
const FieldVector<double,dim>& quadPos = quadPoint.position();
const double integrationElement = element.geometry().integrationElement(quadPos);
// double functionValue = forceTerm(element.geometry().global(quadPos));
// Evaluate all shape function values at this point
// std::vector<FieldVector<double,1> > shapeFunctionValues;
// localFiniteElement.localBasis().evaluateFunction(quadPos, shapeFunctionValues);
const auto jacobian = geometry.jacobianInverseTransposed(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++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
// Actually compute the vector entries
for (size_t p=0; p<nSf; p++)
{
auto localIndex = localView.tree().localIndex(p);
// elementRhs[localIndex] += shapeFunctionValues[p] * forceTerm(quadPos) * quadPoint.weight() * integrationElement;
// auto value = shapeFunctionValues[p]* (sqrt(2.0)*mu(quadPos) *forceTerm(quadPos));
// auto value = -1.0*gradients[p][0] * (sqrt(2.0)*mu(quadPos) *forceTerm(quadPos)); //NEGATIVE!!! TODO
// auto value = -2.0*mu(quadPos)*(sqrt(2.0)*forceTerm(quadPos))*gradients[p][0] ;
// auto value = -1.0*mu(quadPos)*forceTerm(quadPos)*gradients[p][0] ;
// auto value = -2.0*sqrt(2.0)*mu(quadPos)*forceTerm(quadPos)*gradients[p][0]; //#1
auto value = 2.0*sqrt(2.0)*mu(quadPos)*forceTerm(quadPos)*gradients[p][0]; // TEST
// auto value = 2.0*mu(quadPos)*sqrt(2.0)*forceTerm(quadPos)*gradients[p][0];
// auto value = -2.0*mu(quadPos)*sqrt(2.0)*quadPos[1]*gradients[p][0]; //TEST
// std::cout << "value:" << value << std::endl;
// std::cout << "forceTerm(quadPos):" << forceTerm(quadPos) << std::endl;
// std::cout << "quadPos:" << quadPos << std::endl;
// std::cout << "integrationElement:" << integrationElement << std::endl;
// auto value = -1.0*sqrt(2.0)*forceTerm(quadPos)*gradients[p][0] ; //TEST
// auto value = -1.0*mu(quadPos)*sqrt(2.0)*forceTerm(quadPos)*forceTerm(quadPos)*gradients[p][0]; //TEST
elementRhs[localIndex] += value * quadPoint.weight() * integrationElement;
}
}
}
// Get the occupation pattern of the stiffness matrix
template<class Basis>
void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb)
{
nb.resize(basis.size(), basis.size());
auto gridView = basis.gridView();
// A loop over all elements of the grid
auto localView = basis.localView();
for (const auto& element : elements(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,col);
}
}
}
}
/** \brief Assemble the Laplace stiffness matrix on the given grid view */
template<class Basis, class Matrix, class Vector, class LocalFunction, class Force>
void assemblePoissonProblem(const Basis& basis,
Matrix& matrix,
Vector& b,
LocalFunction& muLocal,
const Force& forceTerm,
const double gamma)
{
auto gridView = basis.gridView();
MatrixIndexSet occupationPattern;
getOccupationPattern(basis, occupationPattern);
occupationPattern.exportIdx(matrix);
matrix = 0;
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis.gridView());
auto loadFunctional = localFunction(loadGVF);
b.resize(basis.dimension());
b = 0;
auto localView = basis.localView();
for (const auto& element : elements(gridView))
{
localView.bind(element);
muLocal.bind(element);
loadFunctional.bind(element);
// Dune::Matrix<double> elementMatrix;
Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
// Dune::Matrix<FieldMatrix<double,1,1> > elementMatrix;
assembleElementStiffnessMatrix(localView, elementMatrix, muLocal, gamma);
// std::cout << "elementMatrix.N() "<< elementMatrix.N() << std::endl;
// std::cout << "elementMatrix.M() " << elementMatrix.M() << std::endl;
for(size_t p=0; p<elementMatrix.N(); p++)
{
auto row = localView.index(p);
for (size_t q=0; q<elementMatrix.M(); q++ )
{
auto col = localView.index(q);
matrix[row][col] += elementMatrix[p][q];
}
}
// BlockVector<double> elementRhs;
BlockVector<FieldVector<double,1> > elementRhs;
assembleElementVolumeTerm(localView, elementRhs, muLocal, loadFunctional);
for (size_t p=0; p<elementRhs.size(); p++)
{
// The global index of the p-th vertex of the element
auto row = localView.index(p);
b[row] += elementRhs[p];
}
}
}
template<class Basis, class Vector, class LocalFunc1, class LocalFunc2>
double computeMuGamma(const Basis& basis,
Vector& coeffVector,
const double gamma,
LocalFunc1& mu,
LocalFunc2& Func
)
{
// constexpr int dim = Basis::LocalView::Element::dimension;
double output = 0.0;
double Term1 = 0.0;
double Term2 = 0.0;
double Term11 = 0.0;
constexpr int dim = 2;
// constexpr int dim = 3;
auto localView = basis.localView();
// auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<double>(basis, coeffVector);
// auto localSol = localFunction(solutionFunction);
// auto loadGVF = Dune::Functions::makeGridViewFunction(x3Fun, basis.gridView());
// auto x3Functional = localFunction(loadGVF);
double area = 0.0;
for(const auto& element : elements(basis.gridView()))
{
// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
double elementEnergy1 = 0.0;
double elementEnergy2 = 0.0;
localView.bind(element);
mu.bind(element);
// x3Functional.bind(element);
Func.bind(element);
auto geometry = element.geometry();
const auto& localFiniteElement = localView.tree().finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR); // TODO WARUM HIER KEINE COLOR ?? ERROR
// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
for(const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
// const FieldVector<double,dim>& quadPos = quadPoint.position();
// std::cout << " quadPOS: " << quadPos << std::endl;
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
area += quadPoint.weight() * integrationElement;
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]);
FieldVector<double,dim> tmp(0.0);
// std::cout << "integrationElement :" << integrationElement << std::endl;
// std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
for (size_t i=0; i < nSf; i++)
{
size_t localIdx = localView.tree().localIndex(i); // hier i:leafIdx
size_t globalIdx = localView.index(localIdx);
// printvector(std::cout, gradients[i], "gradients[i]", "--" );
// std::cout << "coeffVector[globalIdx]:" << coeffVector[globalIdx] << std::endl;
tmp[0] += coeffVector[globalIdx]*gradients[i][0];
tmp[1] += coeffVector[globalIdx]*gradients[i][1];
// tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
// tmp[1] += coeffVector[globalIdx]*gradients[i][2];
// printvector(std::cout, tmp, "tmp", "--" );
}
// printvector(std::cout, tmp, "gradient_w", "--" );
// auto value1 = std::pow( ((quadPos[1]/sqrt(2.0))+ (tmp[0]/2.0)) ,2); //TEST
auto q1 = (Func(quadPos)/sqrt(2.0)); //#1
auto q2 = (tmp[0]/2.0); //#1
// auto q2 = (tmp[0]/sqrt(2.0)); // TEST !!!!!!!!!!
// CHECK : BEITRÄGE CHECKEN!!!!
// std::cout << "Beitrag1: " << q2 << std::endl;
// std::cout << "Beitrag2: " << (tmp[1]/(2.0*gamma)) << std::endl;
//
//
//
auto q3 = q1 - q2; // TEST
// auto q3 = q1 + q2; // #1
auto value1 = std::pow(q3,2);
// auto value2 = std::pow( (tmp[1]/(2.0*gamma) ) , 2);
auto value2 = std::pow( (tmp[1]/(sqrt(2.0)*gamma) ) , 2);
// auto value = 2.0*mu(quadPos)*(2.0*value1 + value2);
elementEnergy1 += 2.0*mu(quadPos)* 2.0*value1 * quadPoint.weight() * integrationElement;
elementEnergy2 += 2.0*mu(quadPos)* value2 * quadPoint.weight() * integrationElement;
// std::cout << "output:" << output << std::endl;
}
// std::cout << "elementEnergy:" << elementEnergy << std::endl;
Term1 += elementEnergy1;
Term2 += elementEnergy2;
// std::cout << "output: " << output << std::endl;
}
std::cout << "Term1: " << Term1 << std::endl;
std::cout << "Term2: " << Term2 << std::endl;
output = Term1 + Term2;
std::cout << "output: " << output << std::endl;
std::cout << "Domain-Area: " << area << std::endl;
return (1.0/area) * output;
}
// // -----------------------------------------------------------
/*
template<class Basis, class Vector, class LocalFunc1, class LocalFunc2>
double computeMuGamma(const Basis& basis,
Vector& coeffVector,
const double gamma,
LocalFunc1& mu,
LocalFunc2& Func
)
{
// constexpr int dim = Basis::LocalView::Element::dimension;
double output = 0.0;
constexpr int dim = 2;
// constexpr int dim = 3;
auto localView = basis.localView();
// auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<double>(basis, coeffVector);
// auto localSol = localFunction(solutionFunction);
// auto loadGVF = Dune::Functions::makeGridViewFunction(x3Fun, basis.gridView());
// auto x3Functional = localFunction(loadGVF);
double area = 0.0;
for(const auto& element : elements(basis.gridView()))
{
// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
double elementEnergy = 0.0;
localView.bind(element);
mu.bind(element);
// x3Functional.bind(element);
Func.bind(element);
auto geometry = element.geometry();
const auto& localFiniteElement = localView.tree().finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5;
int orderQR = 2 * (dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 2 * (localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR); // TODO WARUM HIER KEINE COLOR ?? ERROR
// std::cout << " ------------------------- one ELEMENT ------------------------" << std::endl;
for(const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
// const FieldVector<double,dim>& quadPos = quadPoint.position();
// std::cout << " quadPOS: " << quadPos << std::endl;
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
area += quadPoint.weight() * integrationElement;
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]);
// FieldVector<double,dim> tmp = {0.0 , 0.0};
FieldVector<double,dim> tmp(0.0);
// FieldVector<double,dim> tmp = {0.0 ,0.0, 0.0}; //3D-Version
// std::cout << "integrationElement :" << integrationElement << std::endl;
// std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
for (size_t i=0; i < nSf; i++)
{
size_t localIdx = localView.tree().localIndex(i); // hier i:leafIdx
size_t globalIdx = localView.index(localIdx);
// printvector(std::cout, gradients[i], "gradients[i]", "--" );
// std::cout << "coeffVector[globalIdx]:" << coeffVector[globalIdx] << std::endl;
tmp[0] += coeffVector[globalIdx]*gradients[i][0];
tmp[1] += coeffVector[globalIdx]*gradients[i][1];
// tmp[0] += coeffVector[globalIdx]*gradients[i][0]; // 3D-Version
// tmp[1] += coeffVector[globalIdx]*gradients[i][2];
// printvector(std::cout, tmp, "tmp", "--" );
}
// printvector(std::cout, tmp, "gradient_w", "--" );
// auto value1 = std::pow( ((quadPos[1]/sqrt(2.0))+ (tmp[0]/2.0)) ,2); //TEST
auto q1 = (Func(quadPos)/sqrt(2.0));
auto q2 = (tmp[0]/2.0);
// auto q2 = (tmp[0]/sqrt(2.0)); // TEST !!!!!!!!!!
// CHECK : BEITRÄGE CHECKEN!!!!
std::cout << "Beitrag1: " << q2 << std::endl;
std::cout << "Beitrag2: " << (tmp[1]/(2.0*gamma)) << std::endl;
auto q3 = q1 + q2;
auto value1 = std::pow(q3,2);
// auto value1 = std::pow( ((Func(quadPos)/sqrt(2.0)) + (tmp[0]/2.0)) , 2);
auto value2 = std::pow( (tmp[1]/(2.0*gamma) ) , 2);
// auto value2 = std::pow( (tmp[1]/(sqrt(2.0)*gamma) ) , 2); //TEST
// auto value2 = (1.0/(std::pow(gamma,2)))* std::pow(tmp[1],2)/4.0 ; //TEST
auto value = 2.0*mu(quadPos)*(2.0*value1 + value2);
// std::cout << "quadPos[1]:" << quadPos[1]<< std::endl;
// std::cout << "Func(quadPos):" << Func(quadPos) << std::endl;
// std::cout << "sqrt(2.0):" << sqrt(2.0) << std::endl;
// std::cout << "tmp[0]:" << tmp[0] << std::endl;
// std::cout << "tmp[1]:" << tmp[1] << std::endl;
// std::cout << "value1:" << value1 << std::endl;
// std::cout << "value2:" << value2 << std::endl;
// std::cout << "value:" << value << std::endl;
// auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) );
// auto value = 2.0*mu(quadPos)*(2.0* (((x3Functional(quadPos)*x3Functional(quadPos))/2.0) + std::pow( (tmp[0]/2.0) ,2)) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) ); //TEST
// auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) ) + std::pow( (tmp[1]/(2.0*gamma) ) ,2) ; //TEST
// auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2)) + (1.0/gamma)*std::pow( (tmp[1]/2.0) ,2) ; //TEST
// auto value = 2.0*mu(quadPos)*(2.0*std::pow( ((x3Functional(quadPos)/sqrt(2.0))+ (tmp[0]/2.0)) ,2) + (1.0/gamma)*std::pow( (tmp[1]/2.0) ,2) ) ; //TEST
elementEnergy += value * quadPoint.weight() * integrationElement;
// std::cout << "output:" << output << std::endl;
}
// std::cout << "elementEnergy:" << elementEnergy << std::endl;
output += elementEnergy;
// std::cout << "output: " << output << std::endl;
}
std::cout << "Domain-Area: " << area << std::endl;
return (1.0/area) * output;
}
// -------------------------------------------------------------------------
*/
// Check whether two points are equal on R/Z x R/Z
auto equivalent = [](const FieldVector<double,2>& x, const FieldVector<double,2>& 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])) );
};
// // Check whether two points are equal on R/Z x R/Z x R
// auto equivalent3D = [](const FieldVector<double,3>& x, const FieldVector<double,3>& y)
// {
// return ( (FloatCmp::eq(x[0],y[0]) or FloatCmp::eq(x[0]+1,y[0]) or FloatCmp::eq(x[0]-1,y[0]))
// and (FloatCmp::eq(x[1],y[1]) or FloatCmp::eq(x[1]+1,y[1]) or FloatCmp::eq(x[1]-1,y[1]))
// and (FloatCmp::eq(x[2],y[2]))
// );
// };
//
//
int main(int argc, char *argv[])
{
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);
}
//////////////////////////////////
// Generate the grid
//////////////////////////////////
constexpr int dim = 2;
// constexpr int dim = 3; //TEST
// QUAD-GRID
FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0});
// std::array<int, dim> nElements = {128 ,128};
std::array<int, dim> nElements = {64,64};
// std::array<int, dim> nElements = {32 ,32};
// std::array<int, dim> nElements = {16,16};
// std::array<int, dim> nElements = {8,8};
/*
FieldVector<double,dim> lower({-1.0/2.0, -1.0/2.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0/2.0, 1.0/2.0, 1.0/2.0});
std::array<int, dim> nElements = {16,16,16};*/
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
CellGridType grid_CE(lower,upper,nElements);
using GridView = CellGridType::LeafGridView;
const GridView gridView = grid_CE.leafGridView();
using Domain = GridView::Codim<0>::Geometry::GlobalCoordinate;
double gamma = parameterSet.get<double>("gamma", 1.0);
double theta = parameterSet.get<double>("theta", 1.0/4.0);
double mu1 = parameterSet.get<double>("mu1", 10.0);
double beta = parameterSet.get<double>("beta", 2.0);
double mu2 = beta*mu1;
std::cout << "Gamma:" << gamma << std::endl;
std::cout << "Theta:" << theta << std::endl;
std::cout << "mu1:" << mu1 << std::endl;
std::cout << "mu2:" << mu2 << std::endl;
std::cout << "beta:" << beta << std::endl;
auto muTerm = [mu1, mu2, theta] (const Domain& z) {
// if (abs(z[0]) > (theta/2.0))
if (abs(z[0]) >= (theta/2.0))
return mu1;
else
return mu2;
};
// auto materialImp = IsotropicMaterialImp();
// auto muTermT = materialImp.getMu(parameterSet);
auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView);
auto muLocal = localFunction(muGridF);
/////////////////////////////////////////////////////////
// Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
// using Matrix = BCRSMatrix<double>;
// using Vector = BlockVector<double>;
using Vector = BlockVector<FieldVector<double,1> >;
using Matrix = BCRSMatrix<FieldMatrix<double,1,1> >;
Matrix stiffnessMatrix;
Vector b;
/////////////////////////////////////////////////////////
// Assemble the system
/////////////////////////////////////////////////////////
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))
for (const auto& v2 : vertices(gridView))
if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
periodicIndices.unifyIndexPair({gridView.indexSet().index(v1)}, {gridView.indexSet().index(v2)});
auto basis = makeBasis(gridView, Functions::BasisFactory::Experimental::periodic(lagrange<1>(), periodicIndices)); // flatLexicographic()?
// auto x3Func = [](const FieldVector<double,dim>& x){return x[1];}; //2D-Version
// auto x3Func = [](const FieldVector<double,dim>& x){return x[2];}; // 3D-Version
// auto forceTerm = [](const FieldVector<double,dim>& x){return 0.0;}; // gives q3= 2.08333 (close)...
auto forceTerm = [](const FieldVector<double,dim>& x){return x[1];}; //2D-Version
// auto forceTermNEG = [](const FieldVector<double,dim>& x){return -1.0*x[1];}; //2D-Version
// auto forceTerm = [](const FieldVector<double,dim>& x){return x[2];}; // 3D-Version
assemblePoissonProblem(basis, stiffnessMatrix, b, muLocal, forceTerm, gamma);
// printmatrix(std::cout, stiffnessMatrix, "StiffnessMatrix", "--");
// printvector(std::cout, b, "b", "--");
/////////////////////////////////////////////////////////////////////////////
// Write matrix and load vector to files, to be used in later examples
/////////////////////////////////////////////////////////////////////////////
// // { matrix_rhs_writing_begin }
// std::string baseName = "Dune-PeriodicPoisson";
// //+ std::to_string(grid->maxLevel()) + "-refinements";
// storeMatrixMarket(stiffnessMatrix, baseName + "-matrix.mtx");
// storeMatrixMarket(b, baseName + "-rhs.mtx");
// // { matrix_rhs_writing_end }
//
///////////////////////////
// Compute solution
///////////////////////////
// { algebraic_solving_begin }
// Choose an initial iterate that fulfills the Dirichlet conditions
Vector x(basis.size());
x = b;
// Turn the matrix into a linear operator
// MatrixAdapter<Matrix,Vector,Vector> linearOperator(stiffnessMatrix);
//
// // Sequential incomplete LU decomposition as the preconditioner
// SeqILU<Matrix,Vector,Vector> preconditioner(stiffnessMatrix, 1.0); // Relaxation factor
// // Preconditioned conjugate gradient solver
// CGSolver<Vector> cg(linearOperator,
// preconditioner,
// 1e-5, // Desired residual reduction factor
// 50,
// // Maximum number of iterations
// 2);
// // Verbosity of the solver
// // Object storing some statistics about the solving process
// InverseOperatorResult statistics;
std::cout << "------------ CG - Solver ------------" << std::endl;
MatrixAdapter<Matrix, Vector, Vector> op(stiffnessMatrix);
// Sequential incomplete LU decomposition as the preconditioner
SeqILU<Matrix, Vector, Vector> ilu0(stiffnessMatrix,1.0);
int iter = parameterSet.get<double>("iterations_CG", 999);
// Preconditioned conjugate-gradient solver
CGSolver<Vector> solver(op,
ilu0, //NULL,
1e-8, // desired residual reduction factorlack
iter, // maximum number of iterations
2); // verbosity of the solver
InverseOperatorResult statistics;
// Solve!
// cg.apply(x, b, statistics);
solver.apply(x, b, statistics);
// { algebraic_solving_end }
// std::cout << "------------ GMRES - Solver ------------" << std::endl;
// // Turn the matrix into a linear operator
// MatrixAdapter<Matrix, Vector, Vector> op(stiffnessMatrix);
//
// // Fancy (but only) way to not have a preconditioner at all
// Richardson<Vector,Vector> preconditioner(1.0);
//
// // Construct the iterative solver
// RestartedGMResSolver<Vector> solver(
// op, // 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;
//
// // solve for different Correctors (alpha = 1,2,3)
// solver.apply(x, b, statistics);
//
// using SolutionRange = FieldVector<double,dim>;
// using SolutionRange = FieldVector<double,1>;
using SolutionRange = double;
auto solutionFunction = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(basis, x);
// printvector(std::cout, x, "coeffVector", "--" );
std::cout << "Gamma:" << gamma << std::endl;
auto ForceGridF = Dune::Functions::makeGridViewFunction(forceTerm, gridView);
auto ForceLocal = localFunction(ForceGridF);
// double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, x3Func);
// double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, forceTerm);
double mu_gamma = computeMuGamma(basis, x, gamma, muLocal, ForceLocal);
std::cout << "mu_gamma:" << mu_gamma << std::endl;
std::cout.precision(10);
std::cout << "mu_gamma:" << std::fixed << mu_gamma << std::endl;
// TEST "computeMuGamma"
// Vector zeroVec;
// zeroVec.resize(Basis_CE.size());
Vector zeroVec(basis.size());
zeroVec = 0;
std::cout << "TEST computeMuGamma:" << computeMuGamma(basis, zeroVec, gamma, muLocal, ForceLocal)<< std::endl;
std::cout << " --- print analytic solutions --- " << std::endl;
double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
double hm = mu1*(beta/(theta+(1-theta)*beta)) *(1.0/6.0);
double gm = mu1*((1-theta)+theta*beta) *(1.0/6.0);
std::cout << "q1 : " << q1 << std::endl;
std::cout << "q2 : " << q2 << std::endl;
std::cout << "q3 should be between q1 and q2" << std::endl;
// std::cout << "hm : " << hm << std::endl;
// std::cout << "gm : " << gm << std::endl;
if(mu_gamma > q2)
{
std::cout << "mu_gamma > q2!!.. (39) not satisfied" << std::endl;
}
std::cout << "beta : " << beta << std::endl;
std::cout << "beta : " << beta << std::endl;
std::cout << "theta : " << theta << std::endl;
std::cout << "Gamma : " << gamma << std::endl;
std::cout << "mu_gamma:" << mu_gamma << std::endl;
// Output result
VTKWriter<GridView> vtkWriter(gridView);
// vtkWriter.addVertexData(x, "solution"); //--- Anpassen für P2
vtkWriter.addVertexData(
solutionFunction,
VTK::FieldInfo("solution", VTK::FieldInfo::Type::scalar, 1));
// VTK::FieldInfo("solution", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write("Dune-Poisson-Result");
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},{64,64});
using GridViewVTK = VTKGridType::LeafGridView;
const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
std::vector<double> mu_CoeffP0;
Functions::interpolate(scalarP0FeBasis, mu_CoeffP0, muTerm);
auto mu_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, mu_CoeffP0);
VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
MaterialVtkWriter.addCellData(
mu_DGBF_P0,
VTK::FieldInfo("mu_P0", VTK::FieldInfo::Type::scalar, 1));
MaterialVtkWriter.write("MaterialFunctions-level" );
std::cout << "wrote data to file: MaterialFunctions" << std::endl;
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment