#ifndef DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#include <map>
#include <dune/common/parametertree.hh>
// #include <dune/common/float_cmp.hh>
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.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>
#include <dune/microstructure/voigthelper.hh>
using namespace MatrixOperations;
using std::shared_ptr;
using std::make_shared;
using std::fstream;
template <class Basis, class Material> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
class CorrectorComputer {
public:
static const int dimworld = 3; //GridView::dimensionworld;
static const int dim = Basis::GridView::dimension; //const int dim = Domain::dimension;
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using VectorRT = Dune::FieldVector< double, dimworld>;
using MatrixRT = Dune::FieldMatrix< double, dimworld, dimworld>;
using FuncScalar = std::function< double(const Domain&) >;
using FuncVector = std::function< VectorRT(const Domain&) >;
using Func2Tensor = std::function< MatrixRT(const Domain&) >;
using Func2int = std::function< int(const Domain&) >;
using VectorCT = Dune::BlockVector<Dune::FieldVector<double,1> >;
using MatrixCT = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
using ElementMatrixCT = Dune::Matrix<double>;
protected:
//private:
const Basis& basis_;
const Material& material_;
double gamma_;
fstream& log_; // Output-log
const Dune::ParameterTree& parameterSet_;
// const FuncScalar mu_;
// const FuncScalar lambda_;
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
Dune::FieldVector<double,3> m_1_, m_2_, m_3_; // Corrector m_i coefficient vectors
// (assembled) corrector matrices M_i
std::array<MatrixRT, 3 > mContainer;
const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
// ---- Basis for R_sym^{2x2}
// Could really be constexpr static, but then we would need to write the basis here directly in Voigt notation.
// However, I suspect that our Voigt representation may still change in the future.
const std::array<VoigtVector<double,3>,3 > matrixBasis_;
Func2Tensor x3G_1_ = [] (const Domain& x) {
return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
};
Func2Tensor x3G_2_ = [] (const Domain& x) {
return MatrixRT{{0.0, 0.0, 0.0}, {0.0, 1.0*x[2], 0.0}, {0.0, 0.0, 0.0}};
};
Func2Tensor x3G_3_ = [] (const Domain& x) {
return MatrixRT{{0.0, (1.0/sqrt(2.0))*x[2], 0.0}, {(1.0/sqrt(2.0))*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
const std::array<Func2Tensor, 3> x3MatrixBasisContainer_ = {x3G_1_, x3G_2_, x3G_3_};
// --- Offset between basis indices
const int phiOffset_;
public:
///////////////////////////////
// constructor
///////////////////////////////
CorrectorComputer( const Basis& basis,
const Material& material,
double gamma,
std::fstream& log,
const Dune::ParameterTree& parameterSet)
: basis_(basis),
material_(material),
gamma_(gamma),
log_(log),
parameterSet_(parameterSet),
matrixBasis_(std::array<VoigtVector<double,3>,3>{matrixToVoigt(Dune::FieldMatrix<double,3,3>({{1, 0, 0}, {0, 0, 0}, {0, 0, 0}})),
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))}),
phiOffset_(basis.size())
{
// assemble();
//write VTK call here...
}
// -----------------------------------------------------------------
// --- Assemble Corrector problems
void assemble()
{
Dune::Timer StiffnessTimer;
assembleCellStiffness(stiffnessMatrix_);
std::cout << "Stiffness assembly Timer: " << StiffnessTimer.elapsed() << std::endl;
assembleCellLoad(load_alpha1_ ,x3G_1_);
assembleCellLoad(load_alpha2_ ,x3G_2_);
assembleCellLoad(load_alpha3_ ,x3G_3_);
};
///////////////////////////////
// Getter
///////////////////////////////
const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}
Dune::ParameterTree getParameterSet() const {return parameterSet_;}
fstream* getLog(){return &log_;}
double getGamma(){return gamma_;}
shared_ptr<MatrixCT> getStiffnessMatrix(){return make_shared<MatrixCT>(stiffnessMatrix_);}
shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
// shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
// shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
// --- Get Correctors
// shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
// auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
auto getMcontainer(){return mContainer;}
shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}
// shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
auto getMatrixBasiscontainer(){return make_shared<std::array<VoigtVector<double,3>,3 >>(matrixBasis_);}
// auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
// shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
// Get the occupation pattern of the stiffness matrix
void getOccupationPattern(Dune::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)){
Dune::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 < (std::size_t)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
// )
/** \brief Compute the stiffness matrix for one element
*
* \param quadRule The quadrature rule to be used
* \param phaseAtQuadPoint The material phase (i.e., the type of material) at each quadrature point
*/
void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
const Dune::QuadratureRule<double,dim>& quadRule,
const std::vector<int>& phaseAtQuadPoint,
ElementMatrixCT& elementMatrix
)
{
// Local StiffnessMatrix of the form:
// | phi*phi m*phi |
// | phi *m m*m |
auto element = localView.element();
auto geometry = element.geometry();
// using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
elementMatrix = 0;
// auto elasticityTensor = material_.getElasticityTensor();
// LocalBasis-Offset
const int localPhiOffset = localView.size();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
// std::cout << "localView.size(): " << localView.size() << std::endl;
// std::cout << "nSf : " << nSf << std::endl;
int QPcounter= 0;
for (const auto& quadPoint : quadRule)
{
QPcounter++;
const auto& quadPos = quadPoint.position();
const auto geometryJacobianInverse = geometry.jacobianInverse(quadPos);
const auto integrationElement = geometry.integrationElement(quadPos);
std::vector<Dune::FieldMatrix<double,1,dim> > jacobians;
localFiniteElement.localBasis().evaluateJacobian(quadPos, jacobians);
for (size_t i=0; i< jacobians.size(); i++)
jacobians[i] = jacobians[i] * geometryJacobianInverse;
// Compute the scaled deformation gradient for each shape function
std::vector<std::array<VoigtVector<double,dim>,dimworld> > deformationGradient(nSf);
for (size_t i=0; i < nSf; i++)
{
for (size_t k=0; k< dimworld; k++)
{
MatrixRT defGradientV(0);
defGradientV[k] = jacobians[i][0];
deformationGradient[i][k] = symVoigt(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
}
}
// Compute the material phase that the current quadrature point is in
const auto phase = phaseAtQuadPoint[QPcounter-1];
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);
// "phi*phi"-part
for (size_t k=0; k < dimworld; k++)
for (size_t i=0; i < nSf; i++)
{
// auto etmp = material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos));
// auto etmp = elasticityTensor(defGradientU,element.geometry().global(quadPos));
// auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos);
// printmatrix(std::cout, etmp, "etmp", "--");
// double energyDensity= scalarProduct(etmp,defGradientV);
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),defGradientV);
// auto tmp_1 = scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
// if (abs(tmp_1-tmp_2)>1e-2)
// {
// std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
// std::cout << "scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))<< std::endl;
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)<< std::endl;
// }
// Python::Module module = Python::import("material_neukamm");
// auto PythonindicatorFunction = Python::make_function<int>(module.get("f"));
// if (PythonindicatorFunction(element.geometry().global(quadPos)) != material_.applyIndicatorFunction(element.geometry().global(quadPos)))
// {
// std::cout <<" element.geometry().global(quadPos): " << element.geometry().global(quadPos) << std::endl;
// std::cout << "PythonindicatorFunction(element.geometry().global(quadPos)): " << PythonindicatorFunction(element.geometry().global(quadPos)) << std::endl;
// // std::cout << "mu(quadPos):" << mu(quadPos) << std::endl;
// std::cout << "indicatorFunction(element.geometry().global(quadPos)): " << indicatorFunction(element.geometry().global(quadPos)) << std::endl;
// std::cout << "material_.applyIndicatorFunction(element.geometry().global(quadPos)): " << material_.applyIndicatorFunction(element.geometry().global(quadPos)) << std::endl;
// std::cout << "indicatorFunction_material_1: " << indicatorFunction_material_1(element.geometry().global(quadPos)) << std::endl;
// }
// double energyDensity= scalarProduct(material_.ElasticityTensorGlobal(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
// std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
// }
// std::cout << "--------------- \n";
// printmatrix(std::cout, defGradientU, "defGradientU", "--");
// printmatrix(std::cout, material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos))", "--");
// printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
// double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.localElasticityTensor(defGradientU,quadPos,element),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal(defGradientU,quadPos),defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
// double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
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++)
{
// auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) ;
// auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
// if (abs(tmp_1-tmp_2)>1e-2)
// {
// std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
// std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
// }
// double energyDensityGphi = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
double energyDensityGphi = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),deformationGradient[j][l]);
// double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
// std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) <<std::endl;
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(defGradientV));
// double energyDensityGphi= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(defGradientV));
// double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
// double energyDensityGphi= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),defGradientV);
// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
// double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
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;
// auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
// auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
// if (abs(tmp_1-tmp_2)>1e-2)
// {
// std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
// std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])):" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])) << std::endl;
// std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])<< std::endl;
// }
// double energyDensityGG = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
double energyDensityGG = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),matrixBasis_[n]);
// double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
// double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(basisContainer[n]));
// double energyDensityGG= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),basisContainer[n]);
// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
// double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST
elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
// 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
// )
template<class Vector, class LocalForce>
void computeElementLoadVector( const typename Basis::LocalView& localView,
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>;
// auto elasticityTensor = material_.getElasticityTensor();
// auto indicatorFunction = material_.getIndicatorFunction();
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
localIndicatorFunction.bind(element);
// // auto indicatorFunction = material_.getLocalIndicatorFunction();
// auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
// auto indicatorFunction = localFunction(indicatorFunctionGVF);
// indicatorFunction.bind(element);
// Func2int indicatorFunction = *material_.getIndicatorFunction();
// 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();
// 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 = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
// std::cout << "Quad-Rule order used: " << orderQR << std::endl;
for (const auto& quadPoint : quad)
{
const Dune::FieldVector<double,dim>& quadPos = quadPoint.position();
const auto geometryJacobianInverse = geometry.jacobianInverse(quadPos);
const double integrationElement = geometry.integrationElement(quadPos);
std::vector<Dune::FieldMatrix<double,1,dim> > jacobians;
localFiniteElement.localBasis().evaluateJacobian(quadPos, jacobians);
for (size_t i=0; i< jacobians.size(); i++)
jacobians[i] = jacobians[i] * geometryJacobianInverse;
//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 tmpDefGradientV(0);
tmpDefGradientV[k][0] = jacobians[i][0][0]; // Y
tmpDefGradientV[k][1] = jacobians[i][0][1]; // X2
// tmpDefGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2]; //
tmpDefGradientV[k][2] = jacobians[i][0][2]; // X3
VoigtVector<double,3> defGradientV = symVoigt(crossSectionDirectionScaling((1.0/gamma_),tmpDefGradientV));
// double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),defGradientV);
// double energyDensity= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
// double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
// double energyDensity= scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(defGradientV));
// double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),defGradientV);
// double energyDensity= scalarProduct(material_.applyElasticityTensorLocal((-1.0)*forceTerm(quadPos),quadPos),defGradientV);
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV );
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
}
// "f*m"-part
for (size_t m=0; m<3; m++)
{
// double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
double energyDensityfG= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixBasis_[m]);
// double energyDensityfG= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
// double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(basisContainer[m]));
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),basisContainer[m]);
// double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),basisContainer[m]);
// double energyDensityfG = scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),basisContainer[m]);
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
// double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
// std::cout << "energyDensityfG * quadPoint.weight() * integrationElement: " << energyDensityfG * quadPoint.weight() * integrationElement << std::endl;
}
}
}
void assembleCellStiffness(MatrixCT& matrix)
{
std::cout << "assemble Stiffness-Matrix begins." << std::endl;
Dune::MatrixIndexSet occupationPattern;
getOccupationPattern(occupationPattern);
occupationPattern.exportIdx(matrix);
matrix = 0;
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
// A cache for the element stiffness matrices, if desired
bool cacheElementMatrices = true;
std::map<std::vector<int>, ElementMatrixCT> elementMatrixCache;
// 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 auto localPhiOffset = localView.size();
// Dune::Timer Time;
//std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
// Determine a suitable quadrature rule
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quadRule = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
// Determine the material phase at each quadrature point
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
localIndicatorFunction.bind(element);
std::vector<int> phaseAtQuadPoint(quadRule.size());
for (std::size_t i=0; i<quadRule.size(); ++i)
phaseAtQuadPoint[i] = localIndicatorFunction(quadRule[i].position());
ElementMatrixCT elementMatrix;
// We know that all elements of the grid have the same geometry. Therefore, the element stiffness matrices
// depend only on the material phases at the quadrature points. In practice, a few particular phase distributions
// dominate (identical phase at all quadrature points). This makes caching the element matrices worthwhile.
//
// In principle the cache implemented here can get arbitrarily large. If that ever becomes a problem
// we need to implement a smarter cache.
if (cacheElementMatrices)
{
auto cachedMatrix = elementMatrixCache.find(phaseAtQuadPoint);
if (cachedMatrix != elementMatrixCache.end())
{
// This copies the matrix. A smarter implementation would avoid this.
elementMatrix = (*cachedMatrix).second;
}
else
{
computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
// This copies the matrix again.
elementMatrixCache.insert({phaseAtQuadPoint, elementMatrix});
}
}
else // No caching
{
computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
}
// 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:
if(parameterSet_.get<bool>("print_debug", false))
{
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 auto 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, 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();
Dune::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)
Dune::FieldVector<std::size_t,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 = Dune::Functions::makeDiscreteGlobalBasisFunction<Dune::FieldVector<double,dim>>(basis_,coeffVector);
auto localfun = localFunction(GVFunction);
auto localView = basis_.localView();
Dune::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 auto& quad = Dune::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)
{
// Subtract 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;
if(parameterSet_.get<bool>("print_conditionNumber", false))
{
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_;
std::cout << "start corrector solver..." << std::endl;
Dune::Timer SolverTimer;
if (Solvertype==1) // CG - SOLVER
{
std::cout << "------------ CG - Solver ------------" << std::endl;
Dune::MatrixAdapter<MatrixCT, VectorCT, VectorCT> op(stiffnessMatrix_);
// Sequential incomplete LU decomposition as the preconditioner
Dune::SeqILU<MatrixCT, VectorCT, VectorCT> ilu0(stiffnessMatrix_,1.0);
int iter = parameterSet_.get<double>("iterations_CG", 999);
// Preconditioned conjugate-gradient solver
Dune::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
Dune::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
Dune::MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_);
// Fancy (but only) way to not have a preconditioner at all
Dune::Richardson<VectorCT,VectorCT> preconditioner(1.0);
// Construct the iterative solver
Dune::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
Dune::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
Dune::SPQR<MatrixCT> sPQR(stiffnessMatrix_);
sPQR.setVerbosity(1);
Dune::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 << "--- Corrector computation finished ---" << std::endl;
std::cout << "Time required 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) )
for (auto& matrix : mContainer)
matrix = 0;
for(size_t i=0; i<3; i++)
{
mContainer[0] += m_1_[i]*voigtToMatrix(matrixBasis_[i]);
mContainer[1] += m_2_[i]*voigtToMatrix(matrixBasis_[i]);
mContainer[2] += m_3_[i]*voigtToMatrix(matrixBasis_[i]);
}
std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
printmatrix(std::cout, mContainer[0], "Corrector-Matrix M_1", "--");
printmatrix(std::cout, mContainer[1], "Corrector-Matrix M_2", "--");
printmatrix(std::cout, mContainer[2], "Corrector-Matrix M_3", "--");
log_ << "---------- OUTPUT ----------" << std::endl;
log_ << " --------------------" << std::endl;
log_ << "Corrector-Matrix M_1: \n" << mContainer[0] << std::endl;
log_ << " --------------------" << std::endl;
log_ << "Corrector-Matrix M_2: \n" << mContainer[1] << std::endl;
log_ << " --------------------" << std::endl;
log_ << "Corrector-Matrix M_3: \n" << mContainer[2] << 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 << " --- subtracting 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 baseName = parameterSet_.get("baseName", "CellProblem-result");
std::string vtkOutputName = parameterSet_.get("outputPath", "../../outputs") + "/" + baseName;
std::cout << "vtkOutputName:" << vtkOutputName << std::endl;
int subsamplingRefinement = parameterSet_.get<int>("subsamplingRefinement", 2);
Dune::SubsamplingVTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView(), Dune::refinementLevels(subsamplingRefinement));
// Dune::VTKWriter<typename Basis::GridView> vtkWriter(basis_.gridView());
vtkWriter.addVertexData(
Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_1_),
Dune::VTK::FieldInfo("Corrector phi_1 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
vtkWriter.addVertexData(
Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_2_),
Dune::VTK::FieldInfo("Corrector phi_2 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
vtkWriter.addVertexData(
Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_, phi_3_),
Dune::VTK::FieldInfo("Corrector phi_3 level"+ std::to_string(level) , Dune::VTK::FieldInfo::Type::vector, dim));
vtkWriter.write(vtkOutputName + "-level"+ std::to_string(level));
std::cout << "wrote Corrector-VTK data to file: " + vtkOutputName + "-level" + std::to_string(level) << std::endl;
}
// -----------------------------------------------------------------
// --- Compute norms of the corrector functions:
void computeNorms()
{
computeL2Norm();
computeL2SymGrad();
std::cout<< "Frobenius-Norm of M1_: " << mContainer[0].frobenius_norm() << std::endl;
std::cout<< "Frobenius-Norm of M2_: " << mContainer[1].frobenius_norm() << std::endl;
std::cout<< "Frobenius-Norm of M3_: " << mContainer[2].frobenius_norm() << std::endl;
}
void computeL2Norm()
{
// IMPLEMENTATION with makeDiscreteGlobalBasisFunction
double error_1 = 0.0;
double error_2 = 0.0;
double error_3 = 0.0;
auto localView = basis_.localView();
auto GVFunc_1 = Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_);
auto GVFunc_2 = Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_);
auto GVFunc_3 = Dune::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 auto& quad = Dune::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(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
auto GVFunc_2 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
auto GVFunc_3 = derivative(Dune::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();
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );
const auto& quad = Dune::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()
{
Dune::Timer orthogonalityTimer;
std::cout << "Check Orthogonality..." << std::endl;
auto localView = basis_.localView();
auto GVFunc_1 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_1_));
auto GVFunc_2 = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis_,phi_2_));
auto GVFunc_3 = derivative(Dune::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 localIndicatorFunction = material_.getLocalIndicatorFunction();
// auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
// auto mu = localFunction(muGridF);
// auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
// auto lambda= localFunction(lambdaGridF);
for(int a=0; a<3; a++)
for(int b=0; b<3; b++)
{
double energy = 0.0;
auto DerPhi1 = *phiDerContainer[a];
auto DerPhi2 = *phiDerContainer[b];
auto matrixFieldGGVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer_[a], basis_.gridView());
auto matrixFieldG = localFunction(matrixFieldGGVF);
// auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
// auto matrixFieldB = localFunction(matrixFieldBGVF);
// using GridView = typename Basis::GridView;
// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
// using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
// std::cout << "Press key.." << std::endl;
// std::cin.get();
// TEST
// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
for (const auto& e : elements(basis_.gridView()))
{
localView.bind(e);
matrixFieldG.bind(e);
DerPhi1.bind(e);
DerPhi2.bind(e);
// mu.bind(e);
// lambda.bind(e);
localIndicatorFunction.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 auto& quad = Dune::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];
const auto strain1 = symVoigt(crossSectionDirectionScaling(1.0/gamma_, DerPhi1(quadPos)));
auto G = matrixFieldG(quadPos);
// auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
auto tmp = matrixToVoigt(G + mContainer[a]) + strain1;
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),symVoigt(Chi));
// double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
// double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
// elementEnergy += strain1 * quadPoint.weight() * integrationElement;
//elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
//higherPrecEnergy += elementEnergy_HP;
}
if(parameterSet_.get<bool>("print_debug", false))
std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl;
if(energy > 1e-1)
std::cout << "WARNING: check_Orthogonality failed! apparently orthogonality (75) not satisfied on discrete level" << std::endl;
}
std::cout << "Time required for orthogonality check: " << orthogonalityTimer.elapsed() << std::endl;
return 0;
}
// --- Check whether stiffness matrix is symmetric
void checkSymmetry()
{
std::cout << " Check Symmetry of stiffness matrix..." << std::endl;
auto localView = basis_.localView();
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
const auto localPhiOffset = localView.size();
for(size_t i=0; i<localPhiOffset; i++)
for(size_t j=0; j<localPhiOffset; j++ )
{
auto row = localView.index(i);
auto col = localView.index(j);
if(abs( stiffnessMatrix_[row][col] - stiffnessMatrix_[col][row]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
for(size_t i=0; i<localPhiOffset; i++)
for(size_t m=0; m<3; m++)
{
auto row = localView.index(i);
if(abs( stiffnessMatrix_[row][phiOffset_+m] - stiffnessMatrix_[phiOffset_+m][row]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
for(size_t m=0; m<3; m++ )
for(size_t n=0; n<3; n++ )
{
if(abs(stiffnessMatrix_[phiOffset_+m][phiOffset_+n] - stiffnessMatrix_[phiOffset_+n][phiOffset_+m]) > 1e-12 )
std::cout << "STIFFNESS MATRIX NOT SYMMETRIC!!!" << std::endl;
}
}
std::cout << "--- Symmetry test passed ---" << std::endl;
}
}; // end class
#endif