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

Outsource Corrector-Computation to class

parent d0384077
No related branches found
No related tags found
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
dune/microstructure/CorrectorComputer.hh 0 → 100644
+ 1279
− 0
View file @ 60cb7db8
#ifndef DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#include <dune/common/parametertree.hh>
#include <dune/common/float_cmp.hh>
#include <dune/istl/matrixindexset.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
#include <dune/solvers/solvers/umfpacksolver.hh>
using namespace Dune;
// using namespace Functions;
using namespace MatrixOperations;
using std::shared_ptr;
using std::make_shared;
using std::fstream;
template <class Basis> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
class CellAssembler {
public:
static const int dimworld = 3; //GridView::dimensionworld;
static const int dim = Basis::GridView::dimension; //const int dim = Domain::dimension;
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using ScalarRT = FieldVector< double, 1>;
using VectorRT = FieldVector< double, dimworld>;
using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
using FuncScalar = std::function< ScalarRT(const Domain&) >;
using FuncVector = std::function< VectorRT(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using VectorCT = BlockVector<FieldVector<double,1> >;
using MatrixCT = BCRSMatrix<FieldMatrix<double,1,1> >;
using ElementMatrixCT = Matrix<FieldMatrix<double,1,1> >;
using HierarchicVectorView = Dune::Functions::HierarchicVectorWrapper< VectorCT, double>;
protected:
//private:
const Basis& basis_;
fstream& log_; // Output-log
const ParameterTree& parameterSet_;
const FuncScalar mu_;
const FuncScalar lambda_;
double gamma_;
MatrixCT stiffnessMatrix_;
VectorCT load_alpha1_,load_alpha2_,load_alpha3_; //right-hand side(load) vectors
VectorCT x_1_, x_2_, x_3_; // (all) Corrector coefficient vectors
VectorCT phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
FieldVector<double,3> m_1_, m_2_, m_3_; // Corrector m_i coefficient vectors
MatrixRT M1_, M2_, M3_; // (assembled) corrector matrices M_i
const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
// ---- Basis for R_sym^{2x2}
MatrixRT G1_ {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G2_ {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
MatrixRT G3_ {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer_ = {G1_, G2_, G3_};
Func2Tensor x3G_1_ = [] (const Domain& x) {
return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
};
Func2Tensor x3G_2_ = [] (const Domain& x) {
return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
};
Func2Tensor x3G_3_ = [] (const Domain& x) {
return MatrixRT{{0.0, (1.0/sqrt(2.0))*x[2], 0.0}, {(1.0/sqrt(2.0))*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
const std::array<Func2Tensor, 3> x3MatrixBasis_ = {x3G_1_, x3G_2_, x3G_3_};
// --- Offset between basis indices
const int phiOffset_;
public:
///////////////////////////////
// constructor
///////////////////////////////
CellAssembler( const Basis& basis,
const FuncScalar& mu,
const FuncScalar& lambda,
double gamma,
std::fstream& log,
const ParameterTree& parameterSet)
: basis_(basis),
mu_(mu),
lambda_(lambda),
gamma_(gamma),
log_(log),
parameterSet_(parameterSet),
phiOffset_(basis.size())
{
assemble();
// if (parameterSet.get<bool>("stiffnessmatrix_cellproblem_to_csv"))
// csvSystemMatrix();
// if (parameterSet.get<bool>("rhs_cellproblem_to_csv"))
// csvRHSs();
// if (parameterSet.get<bool>("rhs_cellproblem_to_vtk"))
// vtkLoads();
}
// -----------------------------------------------------------------
// --- Assemble Corrector problems
void assemble()
{
assembleCellStiffness(stiffnessMatrix_);
// // lateral stretch strain E_U1 = e1 ot e1 = MatrixRT{{1, 0, 0}, {0, 0, 0}, {0, 0, 0}};
// Func2Tensor neg_E_a = [] (const Domain& z) { return biotStrainApprox({-1, 0, 0}, {0, 0, 0}, {0, z[dim-2], z[dim-1]}); };
// // twist in x2-x3 plane E_K1 = (e1 x (0,x2,x3)^T ot e1) = MatrixRT{{0,-z[2], z[1]}, {0, 0 , 0 }, {0,0,0}};
// Func2Tensor neg_E_K1 = [] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {-1, 0, 0}, {0, z[dim-2], z[dim-1]}); };
// // bend strain in x1-x2 plane E_K2 = (e2 x (0,x2,x3)^T ot e1) = MatrixRT{{z[2], 0, 0}, {0, 0, 0}, {0, 0, 0}};
// Func2Tensor neg_E_K2 = [] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {0, -1, 0}, {0, z[dim-2], z[dim-1]}); };
// // bend strain in x1-x3 plane E_K3 = (e3 x (0,x2,x3)^T ot e1) = MatrixRT{{-z[1], 0, 0}, {0, 0, 0}, {0, 0, 0}};
// Func2Tensor neg_E_K3 = [] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {0, 0, -1}, {0, z[dim-2], z[dim-1]}); };
// assembleLoadVector(neg_E_a, load_a_);
// //log_ << "\n\n\n\n\n\n";
// assembleLoadVector(neg_E_K1, load_K1_);
// assembleLoadVector(neg_E_K2, load_K2_);
// assembleLoadVector(neg_E_K3, load_K3_);
/////////////////////////////////////////////////////////
// Define "rhs"-functions
/////////////////////////////////////////////////////////
// Func2Tensor x3G_1neg = [x3G_1_] (const Domain& x) {return -1.0*x3G_1_(x);};
// Func2Tensor x3G_2neg = [x3G_2_] (const Domain& x) {return -1.0*x3G_2_(x);};
// Func2Tensor x3G_3neg = [x3G_3_] (const Domain& x) {return -1.0*x3G_3_(x);};
// Func2Tensor x3G_1neg = [] (const Domain& x) {return -1.0*x3G_1_(x);};
// Func2Tensor x3G_2neg = [] (const Domain& x) {return -1.0*x3G_2_(x);};
// Func2Tensor x3G_3neg = [] (const Domain& x) {return -1.0*x3G_3_(x);};
// assembleCellLoad(load_alpha1_ ,x3G_1neg);
// assembleCellLoad(load_alpha2_ ,x3G_2neg);
// assembleCellLoad(load_alpha3_ ,x3G_3neg);
assembleCellLoad(load_alpha1_ ,x3G_1_);
assembleCellLoad(load_alpha2_ ,x3G_2_);
assembleCellLoad(load_alpha3_ ,x3G_3_);
};
///////////////////////////////
// getter
///////////////////////////////
const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}
// std::map<int,int> getIndexMapPeriodicBoundary(){return perIndexMap_;}
ParameterTree getParameterSet() const {return parameterSet_;}
fstream* getLog(){return &log_;}
double getGamma(){return gamma_;}
shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
// Get the occupation pattern of the stiffness matrix
void getOccupationPattern(MatrixIndexSet& nb)
{
// OccupationPattern:
// | phi*phi m*phi |
// | phi *m m*m |
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
nb.resize(basis_.size()+3, basis_.size()+3);
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
for (size_t i=0; i<localView.size(); i++)
{
// The global index of the i-th vertex of the element
auto row = localView.index(i);
for (size_t j=0; j<localView.size(); j++ )
{
// The global index of the j-th vertex of the element
auto col = localView.index(j);
nb.add(row[0],col[0]); // nun würde auch nb.add(row,col) gehen..
}
}
for (size_t i=0; i<localView.size(); i++)
{
auto row = localView.index(i);
for (size_t j=0; j<3; j++ )
{
nb.add(row,phiOffset+j); // m*phi part of StiffnessMatrix
nb.add(phiOffset+j,row); // phi*m part of StiffnessMatrix
}
}
for (size_t i=0; i<3; i++ )
for (size_t j=0; j<3; j++ )
{
nb.add(phiOffset+i,phiOffset+j); // m*m part of StiffnessMatrix
}
}
//////////////////////////////////////////////////////////////////
// setOneBaseFunctionToZero
//////////////////////////////////////////////////////////////////
if(parameterSet_.get<bool>("set_oneBasisFunction_Zero ", true)){
FieldVector<int,3> row;
unsigned int arbitraryLeafIndex = parameterSet_.get<unsigned int>("arbitraryLeafIndex", 0);
unsigned int arbitraryElementNumber = parameterSet_.get<unsigned int>("arbitraryElementNumber", 0);
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
const auto nSf = localFiniteElement.localBasis().size();
assert(arbitraryLeafIndex < nSf );
assert(arbitraryElementNumber < basis_.gridView().size(0)); // "arbitraryElementNumber larger than total Number of Elements"
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
for (int k = 0; k<3; k++)
nb.add(row[k],row[k]);
}
std::cout << "finished occupation pattern\n";
}
template<class localFunction1, class localFunction2>
void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
ElementMatrixCT& elementMatrix,
const localFunction1& mu,
const localFunction2& lambda
)
{
// Local StiffnessMatrix of the form:
// | phi*phi m*phi |
// | phi *m m*m |
auto element = localView.element();
auto geometry = element.geometry();
// using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
elementMatrix = 0;
// LocalBasis-Offset
const int localPhiOffset = localView.size();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// std::cout << "localView.size(): " << localView.size() << std::endl;
// std::cout << "nSf : " << nSf << std::endl;
///////////////////////////////////////////////
// Basis for R_sym^{2x2} // wird nicht als Funktion benötigt da konstant...
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}};
MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
// int orderQR = 3;
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
// double elementContribution = 0.0;
// std::cout << "Print QuadratureOrder:" << orderQR << std::endl; //TEST`
int QPcounter= 0;
for (const auto& quadPoint : quad)
{
QPcounter++;
const auto& quadPos = quadPoint.position();
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
const auto integrationElement = geometry.integrationElement(quadPos);
std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
localFiniteElement.localBasis().evaluateJacobian(quadPos,
referenceGradients);
// Compute the shape function gradients on the grid element
std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
for (size_t l=0; l< dimworld; l++)
for (size_t j=0; j < nSf; j++ )
{
size_t row = localView.tree().child(l).localIndex(j);
// (scaled) Deformation gradient of the test basis function
MatrixRT defGradientV(0);
defGradientV[l][0] = gradients[j][0]; // Y
defGradientV[l][1] = gradients[j][1]; //X2
// defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
defGradientV[l][2] = gradients[j][2]; //X3
defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
// "phi*phi"-part
for (size_t k=0; k < dimworld; k++)
for (size_t i=0; i < nSf; i++)
{
// (scaled) Deformation gradient of the ansatz basis function
MatrixRT defGradientU(0);
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
// defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
defGradientU[k][2] = gradients[i][2]; //X3
// printmatrix(std::cout, defGradientU , "defGradientU", "--");
defGradientU = crossSectionDirectionScaling((1.0/gamma_),defGradientU);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
// double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
size_t col = localView.tree().child(k).localIndex(i);
elementMatrix[row][col] += energyDensity * quadPoint.weight() * integrationElement;
}
// "m*phi" & "phi*m" - part
for( size_t m=0; m<3; m++)
{
double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
elementMatrix[localPhiOffset+m][row] += value;
}
}
// "m*m"-part
for(size_t m=0; m<3; m++) //TODO ist symmetric.. reicht die hälfte zu berechnen!!!
for(size_t n=0; n<3; n++)
{
// std::cout << "QPcounter: " << QPcounter << std::endl;
// std::cout << "m= " << m << " n = " << n << std::endl;
// printmatrix(std::cout, basisContainer[m] , "basisContainer[m]", "--");
// printmatrix(std::cout, basisContainer[n] , "basisContainer[n]", "--");
// std::cout << "integrationElement:" << integrationElement << std::endl;
// std::cout << "quadPoint.weight(): "<< quadPoint.weight() << std::endl;
// std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
// std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST
elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
// std::cout << "energyDensityGG:" << energyDensityGG<< std::endl;
// std::cout << "product " << energyDensityGG * quadPoint.weight() * integrationElement << std::endl;
// printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
}
}
// std::cout << "Number of QuadPoints:" << QPcounter << std::endl;
// printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
}
// Compute the source term for a single element
// < L (sym[D_gamma*nabla phi_i] + M_i ), x_3G_alpha >
template<class LocalFunction1, class LocalFunction2, class Vector, class LocalForce>
void computeElementLoadVector( const typename Basis::LocalView& localView,
LocalFunction1& mu,
LocalFunction2& lambda,
Vector& elementRhs,
const LocalForce& forceTerm
)
{
// | f*phi|
// | --- |
// | f*m |
// using Element = typename LocalView::Element;
const auto element = localView.element();
const auto geometry = element.geometry();
// constexpr int dim = Element::dimension;
// constexpr int dimworld = dim;
// using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
// Set of shape functions for a single element
const auto& localFiniteElement= localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
elementRhs.resize(localView.size() +3);
elementRhs = 0;
// LocalBasis-Offset
const int localPhiOffset = localView.size();
///////////////////////////////////////////////
// Basis for R_sym^{2x2}
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}};
MatrixRT G_3 {{0.0, 1.0/sqrt(2.0), 0.0}, {1.0/sqrt(2.0), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
// int orderQR = 3;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
// std::cout << "Quad-Rule order used: " << orderQR << std::endl;
for (const auto& quadPoint : quad)
{
const FieldVector<double,dim>& quadPos = quadPoint.position();
const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
std::vector<FieldMatrix<double,1,dim> > referenceGradients;
localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
std::vector< FieldVector< double, dim> > gradients(referenceGradients.size());
for (size_t i=0; i< gradients.size(); i++)
jacobian.mv(referenceGradients[i][0], gradients[i]);
//TEST
// std::cout << "forceTerm(element.geometry().global(quadPos)):" << std::endl;
// std::cout << forceTerm(element.geometry().global(quadPos)) << std::endl;
// std::cout << "forceTerm(quadPos)" << std::endl;
// std::cout << forceTerm(quadPos) << std::endl;
//
// //TEST QUadrature
// std::cout << "quadPos:" << quadPos << std::endl;
// std::cout << "element.geometry().global(quadPos):" << element.geometry().global(quadPos) << std::endl;
// // //
// // //
// std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
// std::cout << "integrationElement(quadPos):" << integrationElement << std::endl;
// std::cout << "mu(quadPos) :" << mu(quadPos) << std::endl;
// std::cout << "lambda(quadPos) :" << lambda(quadPos) << std::endl;
// "f*phi"-part
for (size_t i=0; i < nSf; i++)
for (size_t k=0; k < dimworld; k++)
{
// Deformation gradient of the ansatz basis function
MatrixRT defGradientV(0);
defGradientV[k][0] = gradients[i][0]; // Y
defGradientV[k][1] = gradients[i][1]; // X2
// defGradientV[k][2] = (1.0/gamma_)*gradients[i][2]; //
defGradientV[k][2] = gradients[i][2]; // X3
defGradientV = crossSectionDirectionScaling((1.0/gamma_),defGradientV);
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV );
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
}
// "f*m"-part
for (size_t m=0; m<3; m++)
{
double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
// std::cout << "energyDensityfG * quadPoint.weight() * integrationElement: " << energyDensityfG * quadPoint.weight() * integrationElement << std::endl;
}
}
}
void assembleCellStiffness(MatrixCT& matrix)
{
std::cout << "assemble Stiffness-Matrix begins." << std::endl;
MatrixIndexSet occupationPattern;
getOccupationPattern(occupationPattern);
occupationPattern.exportIdx(matrix);
matrix = 0;
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
auto muLocal = localFunction(muGridF);
auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
auto lambdaLocal = localFunction(lambdaGridF);
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
muLocal.bind(element);
lambdaLocal.bind(element);
const int localPhiOffset = localView.size();
// Dune::Timer Time;
//std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
Matrix<FieldMatrix<double,1,1> > elementMatrix;
computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
// std::cout << "local assembly time:" << Time.elapsed() << std::endl;
//printmatrix(std::cout, elementMatrix, "ElementMatrix", "--");
//std::cout << "elementMatrix.N() : " << elementMatrix.N() << std::endl;
//std::cout << "elementMatrix.M() : " << elementMatrix.M() << std::endl;
//TEST
//Check Element-Stiffness-Symmetry:
for (size_t i=0; i<localPhiOffset; i++)
for (size_t j=0; j<localPhiOffset; j++ )
{
if(abs(elementMatrix[i][j] - elementMatrix[j][i]) > 1e-12 )
std::cout << "ELEMENT-STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
//////////////////////////////////////////////////////////////////////////////
// GLOBAL STIFFNES ASSEMBLY
//////////////////////////////////////////////////////////////////////////////
for (size_t i=0; i<localPhiOffset; i++)
for (size_t j=0; j<localPhiOffset; j++ )
{
auto row = localView.index(i);
auto col = localView.index(j);
matrix[row][col] += elementMatrix[i][j];
}
for (size_t i=0; i<localPhiOffset; i++)
for(size_t m=0; m<3; m++)
{
auto row = localView.index(i);
matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m];
matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
}
for (size_t m=0; m<3; m++ )
for (size_t n=0; n<3; n++ )
matrix[phiOffset+m][phiOffset+n] += elementMatrix[localPhiOffset+m][localPhiOffset+n];
// printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
}
}
void assembleCellLoad(VectorCT& b,
const Func2Tensor& forceTerm
)
{
// std::cout << "assemble load vector." << std::endl;
b.resize(basis_.size()+3);
b = 0;
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
// Transform G_alpha's to GridViewFunctions/LocalFunctions
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis_.gridView());
auto loadFunctional = localFunction(loadGVF);
auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
auto muLocal = localFunction(muGridF);
auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
auto lambdaLocal = localFunction(lambdaGridF);
// int counter = 1;
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
muLocal.bind(element);
lambdaLocal.bind(element);
loadFunctional.bind(element);
const int localPhiOffset = localView.size();
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
VectorCT elementRhs;
// std::cout << "----------------------------------Element : " << counter << std::endl; //TEST
// counter++;
computeElementLoadVector(localView, muLocal, lambdaLocal, elementRhs, loadFunctional);
// computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
// printvector(std::cout, elementRhs, "elementRhs", "--");
// printvector(std::cout, elementRhs, "elementRhs", "--");
//////////////////////////////////////////////////////////////////////////////
// GLOBAL LOAD ASSEMBLY
//////////////////////////////////////////////////////////////////////////////
for (size_t p=0; p<localPhiOffset; p++)
{
auto row = localView.index(p);
b[row] += elementRhs[p];
}
for (size_t m=0; m<3; m++ )
b[phiOffset+m] += elementRhs[localPhiOffset+m];
//printvector(std::cout, b, "b", "--");
}
}
// -----------------------------------------------------------------
// --- Functions for global integral mean equals zero constraint
auto arbitraryComponentwiseIndices(const int elementNumber,
const int leafIdx
)
{
// (Local Approach -- works for non Lagrangian-Basis too)
// Input : elementNumber & localIdx
// Output : determine all Component-Indices that correspond to that basis-function
auto localView = basis_.localView();
FieldVector<int,3> row;
int elementCounter = 0;
for (const auto& element : elements(basis_.gridView()))
{
if(elementCounter == elementNumber) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
{
localView.bind(element);
for (int k = 0; k < 3; k++)
{
auto rowLocal = localView.tree().child(k).localIndex(leafIdx); //Input: LeafIndex! TODO bräuchte hier (Inverse ) Local-to-Leaf Idx Map
row[k] = localView.index(rowLocal);
// std::cout << "rowLocal:" << rowLocal << std::endl;
// std::cout << "row[k]:" << row[k] << std::endl;
}
}
elementCounter++;
}
return row;
}
void setOneBaseFunctionToZero()
{
std::cout << "Setting one Basis-function to zero" << std::endl;
// constexpr int dim = Basis::LocalView::Element::dimension;
unsigned int arbitraryLeafIndex = parameterSet_.template get<unsigned int>("arbitraryLeafIndex", 0);
unsigned int arbitraryElementNumber = parameterSet_.template get<unsigned int>("arbitraryElementNumber", 0);
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
FieldVector<int,3> row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
for (int k = 0; k<dim; k++)
{
load_alpha1_[row[k]] = 0.0;
load_alpha2_[row[k]] = 0.0;
load_alpha3_[row[k]] = 0.0;
auto cIt = stiffnessMatrix_[row[k]].begin();
auto cEndIt = stiffnessMatrix_[row[k]].end();
for (; cIt!=cEndIt; ++cIt)
*cIt = (cIt.index()==row[k]) ? 1.0 : 0.0;
}
}
auto childToIndexMap(const int k)
{
// Input : child/component
// Output : determine all Indices that belong to that component
auto localView = basis_.localView();
std::vector<int> r = { };
// for (int n: r)
// std::cout << n << ","<< std::endl;
// Determine global indizes for each component k = 1,2,3.. in order to subtract correct (component of) integral Mean
// (global) Indices that correspond to component k = 1,2,3
for(const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
const auto& localFiniteElement = localView.tree().child(k).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
for(size_t j=0; j<nSf; j++)
{
auto Localidx = localView.tree().child(k).localIndex(j); // local indices
r.push_back(localView.index(Localidx)); // global indices
}
}
// Delete duplicate elements
// first remove consecutive (adjacent) duplicates
auto last = std::unique(r.begin(), r.end());
r.erase(last, r.end());
// sort followed by unique, to remove all duplicates
std::sort(r.begin(), r.end());
last = std::unique(r.begin(), r.end());
r.erase(last, r.end());
return r;
}
auto integralMean(VectorCT& coeffVector)
{
auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis_,coeffVector);
auto localfun = localFunction(GVFunction);
auto localView = basis_.localView();
FieldVector<double,3> r = {0.0, 0.0, 0.0};
double area = 0.0;
// Compute Area integral & integral of FE-function
for(const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
localfun.bind(element);
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; //TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
for(const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
const double integrationElement = element.geometry().integrationElement(quadPos);
area += quadPoint.weight() * integrationElement;
r += localfun(quadPos) * quadPoint.weight() * integrationElement;
}
}
// std::cout << "Domain-Area: " << area << std::endl;
return (1.0/area) * r ;
}
auto subtractIntegralMean(VectorCT& coeffVector)
{
// Substract correct Integral mean from each associated component function
auto IM = integralMean(coeffVector);
for(size_t k=0; k<dim; k++)
{
//std::cout << "Integral-Mean: " << IM[k] << std::endl;
auto idx = childToIndexMap(k);
for ( int i : idx)
coeffVector[i] -= IM[k];
}
}
// -----------------------------------------------------------------
// --- Solving the Corrector-problem:
auto solve()
{
bool set_oneBasisFunction_Zero = parameterSet_.get<bool>("set_oneBasisFunction_Zero", false);
bool substract_integralMean = false;
if(parameterSet_.get<bool>("set_IntegralZero", false))
{
set_oneBasisFunction_Zero = true;
substract_integralMean = true;
}
// set one basis-function to zero
if(set_oneBasisFunction_Zero)
setOneBaseFunctionToZero();
//TEST: Compute Condition Number (Can be very expensive !)
const bool verbose = true;
const unsigned int arppp_a_verbosity_level = 2;
const unsigned int pia_verbosity_level = 1;
MatrixInfo<MatrixCT> matrixInfo(stiffnessMatrix_,verbose,arppp_a_verbosity_level,pia_verbosity_level);
std::cout << "Get condition number of Stiffness_CE: " << matrixInfo.getCond2(true) << std::endl;
///////////////////////////////////
// --- Choose Solver ---
// 1 : CG-Solver
// 2 : GMRES
// 3 : QR (default)
// 4 : UMFPACK
///////////////////////////////////
unsigned int Solvertype = parameterSet_.get<unsigned int>("Solvertype", 3);
unsigned int Solver_verbosity = parameterSet_.get<unsigned int>("Solver_verbosity", 2);
// --- set initial values for solver
x_1_ = load_alpha1_;
x_2_ = load_alpha2_;
x_3_ = load_alpha3_;
Dune::Timer SolverTimer;
if (Solvertype==1) // CG - SOLVER
{
std::cout << "------------ CG - Solver ------------" << std::endl;
MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);
// Sequential incomplete LU decomposition as the preconditioner
SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_,1.0);
int iter = parameterSet_.get<double>("iterations_CG", 999);
// Preconditioned conjugate-gradient solver
CGSolver<VectorCT> solver(op,
ilu0, //NULL,
1e-8, // desired residual reduction factorlack
iter, // maximum number of iterations
Solver_verbosity,
true // Try to estimate condition_number
); // verbosity of the solver
InverseOperatorResult statistics;
std::cout << "solve linear system for x_1.\n";
solver.apply(x_1_, load_alpha1_, statistics);
std::cout << "solve linear system for x_2.\n";
solver.apply(x_2_, load_alpha2_, statistics);
std::cout << "solve linear system for x_3.\n";
solver.apply(x_3_, load_alpha3_, statistics);
log_ << "Solver-type used: " <<" CG-Solver" << std::endl;
std::cout << "statistics.converged " << statistics.converged << std::endl;
std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
std::cout << "statistics.iterations: " << statistics.iterations << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
else if (Solvertype==2) // GMRES - SOLVER
{
std::cout << "------------ GMRES - Solver ------------" << std::endl;
// Turn the matrix into a linear operator
MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_);
// Fancy (but only) way to not have a preconditioner at all
Richardson<VectorCT,VectorCT> preconditioner(1.0);
// Construct the iterative solver
RestartedGMResSolver<VectorCT> solver(
stiffnessOperator, // Operator to invert
preconditioner, // Preconditioner
1e-10, // Desired residual reduction factor
500, // Number of iterations between restarts,
// here: no restarting
500, // Maximum number of iterations
Solver_verbosity); // Verbosity of the solver
// Object storing some statistics about the solving process
InverseOperatorResult statistics;
// solve for different Correctors (alpha = 1,2,3)
solver.apply(x_1_, load_alpha1_, statistics); //load_alpha1 now contains the corresponding residual!!
solver.apply(x_2_, load_alpha2_, statistics);
solver.apply(x_3_, load_alpha3_, statistics);
log_ << "Solver-type used: " <<" GMRES-Solver" << std::endl;
std::cout << "statistics.converged " << statistics.converged << std::endl;
std::cout << "statistics.condition_estimate: " << statistics.condition_estimate << std::endl;
std::cout << "statistics.iterations: " << statistics.iterations << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
else if ( Solvertype==3)// QR - SOLVER
{
std::cout << "------------ QR - Solver ------------" << std::endl;
log_ << "solveLinearSystems: We use QR solver.\n";
//TODO install suitesparse
SPQR<MatrixCT> sPQR(stiffnessMatrix_);
sPQR.setVerbosity(1);
InverseOperatorResult statisticsQR;
sPQR.apply(x_1_, load_alpha1_, statisticsQR);
std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
sPQR.apply(x_2_, load_alpha2_, statisticsQR);
std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
sPQR.apply(x_3_, load_alpha3_, statisticsQR);
std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
log_ << "Solver-type used: " <<" QR-Solver" << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
else if (Solvertype==4)// UMFPACK - SOLVER
{
std::cout << "------------ UMFPACK - Solver ------------" << std::endl;
log_ << "solveLinearSystems: We use UMFPACK solver.\n";
Dune::Solvers::UMFPackSolver<MatrixCT,VectorCT> solver;
solver.setProblem(stiffnessMatrix_,x_1_,load_alpha1_);
// solver.preprocess();
solver.solve();
solver.setProblem(stiffnessMatrix_,x_2_,load_alpha2_);
// solver.preprocess();
solver.solve();
solver.setProblem(stiffnessMatrix_,x_3_,load_alpha3_);
// solver.preprocess();
solver.solve();
// sPQR.apply(x_1, load_alpha1, statisticsQR);
// std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
// std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
// std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
// sPQR.apply(x_2, load_alpha2, statisticsQR);
// std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
// std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
// std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
// sPQR.apply(x_3, load_alpha3, statisticsQR);
// std::cout << "statistics.converged " << statisticsQR.converged << std::endl;
// std::cout << "statistics.condition_estimate: " << statisticsQR.condition_estimate << std::endl;
// std::cout << "statistics.iterations: " << statisticsQR.iterations << std::endl;
log_ << "Solver-type used: " <<" UMFPACK-Solver" << std::endl;
}
std::cout << "Finished solving Corrector problems!" << std::endl;
std::cout << "Time for solving:" << SolverTimer.elapsed() << std::endl;
////////////////////////////////////////////////////////////////////////////////////
// Extract phi_alpha & M_alpha coefficients
////////////////////////////////////////////////////////////////////////////////////
phi_1_.resize(basis_.size());
phi_1_ = 0;
phi_2_.resize(basis_.size());
phi_2_ = 0;
phi_3_.resize(basis_.size());
phi_3_ = 0;
for(size_t i=0; i<basis_.size(); i++)
{
phi_1_[i] = x_1_[i];
phi_2_[i] = x_2_[i];
phi_3_[i] = x_3_[i];
}
for(size_t i=0; i<3; i++)
{
m_1_[i] = x_1_[phiOffset_+i];
m_2_[i] = x_2_[phiOffset_+i];
m_3_[i] = x_3_[phiOffset_+i];
}
// assemble M_alpha's (actually iota(M_alpha) )
// MatrixRT M_1(0), M_2(0), M_3(0);
M1_ = 0;
M2_ = 0;
M3_ = 0;
for(size_t i=0; i<3; i++)
{
M1_ += m_1_[i]*basisContainer_[i];
M2_ += m_2_[i]*basisContainer_[i];
M3_ += m_3_[i]*basisContainer_[i];
}
std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
printmatrix(std::cout, M1_, "Corrector-Matrix M_1", "--");
printmatrix(std::cout, M2_, "Corrector-Matrix M_2", "--");
printmatrix(std::cout, M3_, "Corrector-Matrix M_3", "--");
log_ << "---------- OUTPUT ----------" << std::endl;
log_ << " --------------------" << std::endl;
log_ << "Corrector-Matrix M_1: \n" << M1_ << std::endl;
log_ << " --------------------" << std::endl;
log_ << "Corrector-Matrix M_2: \n" << M2_ << std::endl;
log_ << " --------------------" << std::endl;
log_ << "Corrector-Matrix M_3: \n" << M3_ << std::endl;
log_ << " --------------------" << std::endl;
if(parameterSet_.get<bool>("write_IntegralMean", false))
{
std::cout << "check integralMean phi_1: " << std::endl;
auto A = integralMean(phi_1_);
for(size_t i=0; i<3; i++)
{
std::cout << "Integral-Mean phi_1 : " << A[i] << std::endl;
}
}
if(substract_integralMean)
{
std::cout << " --- substracting integralMean --- " << std::endl;
subtractIntegralMean(phi_1_);
subtractIntegralMean(phi_2_);
subtractIntegralMean(phi_3_);
subtractIntegralMean(x_1_);
subtractIntegralMean(x_2_);
subtractIntegralMean(x_3_);
//////////////////////////////////////////
// Check Integral-mean again:
//////////////////////////////////////////
if(parameterSet_.get<bool>("write_IntegralMean", false))
{
auto A = integralMean(phi_1_);
for(size_t i=0; i<3; i++)
{
std::cout << "Integral-Mean phi_1 (Check again) : " << A[i] << std::endl;
}
}
}
/////////////////////////////////////////////////////////
// Write Solution (Corrector Coefficients) in Logs
/////////////////////////////////////////////////////////
if(parameterSet_.get<bool>("write_corrector_phi1", false))
{
log_ << "\nSolution of Corrector problems:\n";
log_ << "\n Corrector_phi1:\n";
log_ << x_1_ << std::endl;
}
if(parameterSet_.get<bool>("write_corrector_phi2", false))
{
log_ << "-----------------------------------------------------";
log_ << "\n Corrector_phi2:\n";
log_ << x_2_ << std::endl;
}
if(parameterSet_.get<bool>("write_corrector_phi3", false))
{
log_ << "-----------------------------------------------------";
log_ << "\n Corrector_phi3:\n";
log_ << x_3_ << std::endl;
}
}
// -----------------------------------------------------------------
// --- Write Correctos to VTK:
void writeCorrectorsVTK(const int level)
{
std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/CellProblem-result";
std::cout << "vtkOutputName:" << vtkOutputName << std::endl;
VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
vtkWriter.addVertexData(
Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_1_),
VTK::FieldInfo("Corrector phi_1 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
vtkWriter.addVertexData(
Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_2_),
VTK::FieldInfo("Corrector phi_2 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
vtkWriter.addVertexData(
Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_3_),
VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(vtkOutputName + "-level"+ std::to_string(level));
std::cout << "wrote Corrector-VTK data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;
}
// -----------------------------------------------------------------
// --- Compute norms of the corrector functions:
void computeNorms()
{
computeL2Norm();
computeL2SymGrad();
}
void computeL2Norm()
{
// IMPLEMENTATION with makeDiscreteGlobalBasisFunction
double error_1 = 0.0;
double error_2 = 0.0;
double error_3 = 0.0;
auto localView = basis_.localView();
auto GVFunc_1 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_);
auto GVFunc_2 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_);
auto GVFunc_3 = Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_);
auto localfun_1 = localFunction(GVFunc_1);
auto localfun_2 = localFunction(GVFunc_2);
auto localfun_3 = localFunction(GVFunc_3);
for(const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
localfun_1.bind(element);
localfun_2.bind(element);
localfun_3.bind(element);
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
for(const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
const double integrationElement = element.geometry().integrationElement(quadPos);
error_1 += localfun_1(quadPos)*localfun_1(quadPos) * quadPoint.weight() * integrationElement;
error_2 += localfun_2(quadPos)*localfun_2(quadPos) * quadPoint.weight() * integrationElement;
error_3 += localfun_3(quadPos)*localfun_3(quadPos) * quadPoint.weight() * integrationElement;
}
}
std::cout << "L2-Norm(Corrector phi_1): " << sqrt(error_1) << std::endl;
std::cout << "L2-Norm(Corrector phi_2): " << sqrt(error_2) << std::endl;
std::cout << "L2-Norm(Corrector phi_3): " << sqrt(error_3) << std::endl;
return;
}
void computeL2SymGrad()
{
double error_1 = 0.0;
double error_2 = 0.0;
double error_3 = 0.0;
auto localView = basis_.localView();
auto GVFunc_1 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
auto GVFunc_2 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
auto GVFunc_3 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
auto localfun_1 = localFunction(GVFunc_1);
auto localfun_2 = localFunction(GVFunc_2);
auto localfun_3 = localFunction(GVFunc_3);
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
localfun_1.bind(element);
localfun_2.bind(element);
localfun_3.bind(element);
auto geometry = element.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
const auto integrationElement = geometry.integrationElement(quadPos);
auto scaledSymGrad1 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_1(quadPos)));
auto scaledSymGrad2 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_2(quadPos)));
auto scaledSymGrad3 = sym(crossSectionDirectionScaling(1.0/gamma_, localfun_3(quadPos)));
error_1 += scalarProduct(scaledSymGrad1,scaledSymGrad1) * quadPoint.weight() * integrationElement;
error_2 += scalarProduct(scaledSymGrad2,scaledSymGrad2) * quadPoint.weight() * integrationElement;
error_3 += scalarProduct(scaledSymGrad3,scaledSymGrad3) * quadPoint.weight() * integrationElement;
}
}
std::cout << "L2-Norm(Symmetric scaled gradient of Corrector phi_1): " << sqrt(error_1) << std::endl;
std::cout << "L2-Norm(Symmetric scaled gradient of Corrector phi_2): " << sqrt(error_2) << std::endl;
std::cout << "L2-Norm(Symmetric scaled gradient of Corrector phi_3): " << sqrt(error_3) << std::endl;
return;
}
// -----------------------------------------------------------------
// --- Check Orthogonality relation Paper (75)
auto check_Orthogonality()
{
std::cout << "CHECK ORTHOGONALITY" << std::endl;
auto localView = basis_.localView();
auto GVFunc_1 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
auto GVFunc_2 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
auto GVFunc_3 = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_3_));
auto localfun_1 = localFunction(GVFunc_1);
auto localfun_2 = localFunction(GVFunc_2);
auto localfun_3 = localFunction(GVFunc_3);
const std::array<decltype(localfun_1)*,3> phiDerContainer = {&localfun_1 , &localfun_2 , &localfun_3 };
auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
auto mu = localFunction(muGridF);
auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
auto lambda= localFunction(lambdaGridF);
for(int a=0; a<3; a++)
for(int b=0; b<3; b++)
{
double energy = 0.0;
auto DerPhi1 = *phiDerContainer[a];
auto DerPhi2 = *phiDerContainer[b];
auto matrixFieldGGVF = Dune::Functions::makeGridViewFunction(x3MatrixBasis_[a], basis_.gridView());
auto matrixFieldG = localFunction(matrixFieldGGVF);
// auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
// auto matrixFieldB = localFunction(matrixFieldBGVF);
// using GridView = typename Basis::GridView;
// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
// using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
// std::cout << "Press key.." << std::endl;
// std::cin.get();
// TEST
// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
for (const auto& e : elements(basis_.gridView()))
{
localView.bind(e);
matrixFieldG.bind(e);
DerPhi1.bind(e);
DerPhi2.bind(e);
mu.bind(e);
lambda.bind(e);
double elementEnergy = 0.0;
//double elementEnergy_HP = 0.0;
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
// int orderQR = 0;
// int orderQR = 1;
// int orderQR = 2;
// int orderQR = 3;
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
const double integrationElement = geometry.integrationElement(quadPos);
auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + *mContainer[b];
auto strain1 = DerPhi1(quadPos);
// printmatrix(std::cout, strain1 , "strain1", "--");
//cale with GAMMA
strain1 = crossSectionDirectionScaling(1.0/gamma_, strain1);
strain1 = sym(strain1);
// ADD M
// auto test = strain1 + *M ;
// std::cout << "test:" << test << std::endl;
// for (size_t m=0; m<3; m++ )
// for (size_t n=0; n<3; n++ )
// strain1[m][n] += M[m][n];
auto G = matrixFieldG(quadPos);
// auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
auto tmp = G + *mContainer[a] + strain1;
double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
// elementEnergy += strain1 * quadPoint.weight() * integrationElement;
//elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
//higherPrecEnergy += elementEnergy_HP;
}
std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl;
// TEST
// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
}
return 0;
}
}; // end class
#endif
\ No newline at end of file
src/Cell-Problem-New.cc 0 → 100644
+ 2746
− 0
View file @ 60cb7db8
This diff is collapsed.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment