diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index a2fef077f1c826abb1341e197e3d32fd19530780..d5b7c9b7365b91ea82567d57d50301c23b505e04 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -1,7 +1,6 @@
#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>
@@ -23,7 +22,7 @@ using std::fstream;
template <class Basis> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
-class CellAssembler {
+class CorrectorComputer {
public:
static const int dimworld = 3; //GridView::dimensionworld;
@@ -54,1221 +53,1228 @@ protected:
const ParameterTree& parameterSet_;
const FuncScalar mu_;
- const FuncScalar lambda_;
- double gamma_;
+ const FuncScalar lambda_;
+ double gamma_;
- MatrixCT stiffnessMatrix_;
+ 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
+ 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_};
+ MatrixRT M1_, M2_, M3_; // (assembled) corrector matrices M_i
+ const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
+ const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
- // ---- 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_};
+ // ---- 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) {
+ 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) {
+ 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) {
+ 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_};
+ const std::array<Func2Tensor, 3> x3MatrixBasis_ = {x3G_1_, x3G_2_, x3G_3_};
// --- Offset between basis indices
- const int phiOffset_;
+ 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())
- {
+ ///////////////////////////////
+ // constructor
+ ///////////////////////////////
+ CorrectorComputer( 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();
+ 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();
- }
+ // 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_);
+ // -----------------------------------------------------------------
+ // --- 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]}); };
+ // // 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]}); };
+ // // 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]}); };
+ // // 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
+ /////////////////////////////////////////////////////////
- // 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);
-// 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_);
+ };
- 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_);}
- ///////////////////////////////
- // getter
- ///////////////////////////////
- const shared_ptr<Basis> getBasis() {return make_shared<Basis>(basis_);}
+ // std::map<int,int> getIndexMapPeriodicBoundary(){return perIndexMap_;}
- // std::map<int,int> getIndexMapPeriodicBoundary(){return perIndexMap_;}
+ ParameterTree getParameterSet() const {return parameterSet_;}
- ParameterTree getParameterSet() const {return parameterSet_;}
+ fstream* getLog(){return &log_;}
- fstream* getLog(){return &log_;}
+ double getGamma(){return gamma_;}
- 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<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<VectorCT> getMu(){return make_shared<FuncScalar>(mu_);}
+ shared_ptr<VectorCT> getLambda(){return make_shared<FuncScalar>(lambda_);}
-// 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();
+ // --- Get Correctors
+ shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
+ shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}
+ // shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
+ // shared_ptr<VectorCT> getCorr_K1(){return make_shared<VectorCT>(x_K1_);}
+ // shared_ptr<VectorCT> getCorr_K2(){return make_shared<VectorCT>(x_K2_);}
+ // shared_ptr<VectorCT> getCorr_K3(){return make_shared<VectorCT>(x_K3_);}
- nb.resize(basis_.size()+3, basis_.size()+3);
- for (const auto& element : elements(basis_.gridView()))
+ // Get the occupation pattern of the stiffness matrix
+ void getOccupationPattern(MatrixIndexSet& nb)
{
- localView.bind(element);
- for (size_t i=0; i<localView.size(); i++)
+ // 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()))
{
- // 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++ )
+ 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++)
{
- // 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..
+ 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]);
}
- for (size_t i=0; i<localView.size(); i++)
+ 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)
{
- auto row = localView.index(i);
+ QPcounter++;
+ const auto& quadPos = quadPoint.position();
+ const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
+ const auto integrationElement = geometry.integrationElement(quadPos);
- for (size_t j=0; j<3; j++ )
+ 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++ )
{
- nb.add(row,phiOffset+j); // m*phi part of StiffnessMatrix
- nb.add(phiOffset+j,row); // phi*m part of StiffnessMatrix
+ 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", "--");
+ }
}
- 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
- }
+ // std::cout << "Number of QuadPoints:" << QPcounter << std::endl;
+ // printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
}
- //////////////////////////////////////////////////////////////////
- // 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)
+ // 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
+ )
{
- 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++ )
+ // | 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)
{
- 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
+ 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++)
- 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;
- }
+ // 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
- }
- // "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++)
- {
+ 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
-// 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", "--");
+ 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;
}
+ }
}
-// 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)
+
+
+ void assembleCellStiffness(MatrixCT& matrix)
{
- const FieldVector<double,dim>& quadPos = quadPoint.position();
- const auto jacobian = geometry.jacobianInverseTransposed(quadPos);
- const double integrationElement = geometry.integrationElement(quadPos);
+ std::cout << "assemble Stiffness-Matrix begins." << std::endl;
- std::vector<FieldMatrix<double,1,dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPos, referenceGradients);
+ MatrixIndexSet occupationPattern;
+ getOccupationPattern(occupationPattern);
+ occupationPattern.exportIdx(matrix);
+ matrix = 0;
- 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;
-
+ auto localView = basis_.localView();
+ const int phiOffset = basis_.dimension();
- // "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);
+ auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ auto muLocal = localFunction(muGridF);
+ auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ auto lambdaLocal = localFunction(lambdaGridF);
- 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;
+ 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];
- // "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;
+ // printmatrix(std::cout, matrix, "StiffnessMatrix", "--");
}
}
-}
-void assembleCellStiffness(MatrixCT& matrix)
-{
- std::cout << "assemble Stiffness-Matrix begins." << std::endl;
+ void assembleCellLoad(VectorCT& b,
+ const Func2Tensor& forceTerm
+ )
+ {
+ // std::cout << "assemble load vector." << std::endl;
+ b.resize(basis_.size()+3);
+ b = 0;
- MatrixIndexSet occupationPattern;
- getOccupationPattern(occupationPattern);
- occupationPattern.exportIdx(matrix);
- matrix = 0;
+ auto localView = basis_.localView();
+ const int phiOffset = basis_.dimension();
- 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);
+ 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++)
+
+ // int counter = 1;
+ for (const auto& element : elements(basis_.gridView()))
{
- auto row = localView.index(i);
- matrix[row][phiOffset+m] += elementMatrix[i][localPhiOffset+m];
- matrix[phiOffset+m][row] += elementMatrix[localPhiOffset+m][i];
+ 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", "--");
}
- 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();
+ // -----------------------------------------------------------------
+ // --- 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();
- // Transform G_alpha's to GridViewFunctions/LocalFunctions
- auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis_.gridView());
- auto loadFunctional = localFunction(loadGVF);
+ FieldVector<int,3> row;
+ int elementCounter = 0;
- 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()))
+ {
+ 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;
+ }
-// int counter = 1;
- for (const auto& element : elements(basis_.gridView()))
+ void setOneBaseFunctionToZero()
{
- localView.bind(element);
- muLocal.bind(element);
- lambdaLocal.bind(element);
- loadFunctional.bind(element);
+ std::cout << "Setting one Basis-function to zero" << std::endl;
- 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++)
+ // 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++)
{
- auto row = localView.index(p);
- b[row] += elementRhs[p];
+ 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;
}
- 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()))
+
+
+ auto childToIndexMap(const int k)
{
- if(elementCounter == elementNumber) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
+ // 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 (int k = 0; k < 3; k++)
+ for(size_t j=0; j<nSf; j++)
{
- 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;
+ auto Localidx = localView.tree().child(k).localIndex(j); // local indices
+ r.push_back(localView.index(Localidx)); // global indices
}
}
- elementCounter++;
+ // 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;
}
- 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++)
+ auto integralMean(VectorCT& coeffVector)
{
- 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 GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis_,coeffVector);
+ auto localfun = localFunction(GVFunction);
-auto childToIndexMap(const int k)
-{
- // Input : child/component
- // Output : determine all Indices that belong to that component
- auto localView = basis_.localView();
+ 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();
+ FieldVector<double,3> r = {0.0, 0.0, 0.0};
+ double area = 0.0;
- for(size_t j=0; j<nSf; j++)
+ // Compute Area integral & integral of FE-function
+ for(const auto& element : elements(basis_.gridView()))
{
- auto Localidx = localView.tree().child(k).localIndex(j); // local indices
- r.push_back(localView.index(Localidx)); // global indices
+ 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 ;
}
- // 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()))
+
+
+ auto subtractIntegralMean(VectorCT& coeffVector)
{
- 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);
+ // Substract correct Integral mean from each associated component function
+ auto IM = integralMean(coeffVector);
- for(const auto& quadPoint : quad)
+ for(size_t k=0; k<dim; k++)
{
- const auto& quadPos = quadPoint.position();
- const double integrationElement = element.geometry().integrationElement(quadPos);
- area += quadPoint.weight() * integrationElement;
-
- r += localfun(quadPos) * quadPoint.weight() * integrationElement;
+ //std::cout << "Integral-Mean: " << IM[k] << std::endl;
+ auto idx = childToIndexMap(k);
+ for ( int i : idx)
+ coeffVector[i] -= IM[k];
}
}
- // 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++)
+ // -----------------------------------------------------------------
+ // --- Solving the Corrector-problem:
+ auto solve()
{
- //std::cout << "Integral-Mean: " << IM[k] << std::endl;
- auto idx = childToIndexMap(k);
- for ( int i : idx)
- coeffVector[i] -= IM[k];
- }
-}
+ 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;
-// -----------------------------------------------------------------
-// --- Solving the Corrector-problem:
-auto solve()
-{
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+ 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) )
- 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;
+ // MatrixRT M_1(0), M_2(0), M_3(0);
- }
- ////////////////////////////////////////////////////////////////////////////////////
- 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) )
+ M1_ = 0;
+ M2_ = 0;
+ M3_ = 0;
- // MatrixRT M_1(0), M_2(0), M_3(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];
+ }
- M1_ = 0;
- M2_ = 0;
- M3_ = 0;
+ 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;
+ }
- 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()))
+ // -----------------------------------------------------------------
+ // --- Write Correctos to VTK:
+ void writeCorrectorsVTK(const int level)
{
- localView.bind(element);
- localfun_1.bind(element);
- localfun_2.bind(element);
- localfun_3.bind(element);
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ 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;
+ }
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
+ // -----------------------------------------------------------------
+ // --- Compute norms of the corrector functions:
+ void computeNorms()
+ {
+ computeL2Norm();
+ computeL2SymGrad();
- 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<< "Frobenius-Norm of M1_: " << M1_.frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M2_: " << M2_.frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M3_: " << M3_.frobenius_norm() << std::endl;
}
- 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 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();
-void computeL2SymGrad()
-{
- double error_1 = 0.0;
- double error_2 = 0.0;
- double error_3 = 0.0;
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
- auto localView = basis_.localView();
+ 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;
- 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);
+ return;
+ }
- for (const auto& element : elements(basis_.gridView()))
+ void computeL2SymGrad()
{
- localView.bind(element);
- localfun_1.bind(element);
- localfun_2.bind(element);
- localfun_3.bind(element);
- auto geometry = element.geometry();
+ double error_1 = 0.0;
+ double error_2 = 0.0;
+ double error_3 = 0.0;
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- const auto nSf = localFiniteElement.localBasis().size();
+ auto localView = basis_.localView();
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );
- const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
+ 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& quadPoint : quad)
+ for (const auto& element : elements(basis_.gridView()))
{
- const auto& quadPos = quadPoint.position();
- const auto integrationElement = geometry.integrationElement(quadPos);
+ localView.bind(element);
+ localfun_1.bind(element);
+ localfun_2.bind(element);
+ localfun_3.bind(element);
+ auto geometry = element.geometry();
- 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)));
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ const auto nSf = localFiniteElement.localBasis().size();
- error_1 += scalarProduct(scaledSymGrad1,scaledSymGrad1) * quadPoint.weight() * integrationElement;
- error_2 += scalarProduct(scaledSymGrad2,scaledSymGrad2) * quadPoint.weight() * integrationElement;
- error_3 += scalarProduct(scaledSymGrad3,scaledSymGrad3) * quadPoint.weight() * integrationElement;
+ 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;
}
- 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)
+ // -----------------------------------------------------------------
+ // --- 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 check_Orthogonality()
-{
- std::cout << "CHECK ORTHOGONALITY" << std::endl;
- auto localView = basis_.localView();
+ auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
+ auto mu = localFunction(muGridF);
+ auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
+ auto lambda= localFunction(lambdaGridF);
- 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 };
+ 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) ", "--");
- auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- auto mu = localFunction(muGridF);
- auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- auto lambda= localFunction(lambdaGridF);
+ 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];
- for(int a=0; a<3; a++)
- for(int b=0; b<3; b++)
- {
- double energy = 0.0;
+ auto strain1 = DerPhi1(quadPos);
- 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);
-
+ strain1 = crossSectionDirectionScaling(1.0/gamma_, strain1);
+ strain1 = sym(strain1);
+
- auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + *mContainer[b];
+ auto G = matrixFieldG(quadPos);
+ // auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
+
+ auto tmp = G + *mContainer[a] + strain1;
- 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;
+ }
+ if(parameterSet_.get<bool>("print_debug", false))
+ std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl;
- double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
+ if(energy > 1e-1)
+ std::cout << "WARNING: check Orthogonality! apparently (75) not satisfied on discrete level" << std::endl;
- 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;
}
- return 0;
-}
-
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
new file mode 100644
index 0000000000000000000000000000000000000000..12151857570b04bc46aeec54cec2a697a8e2387e
--- /dev/null
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -0,0 +1,254 @@
+#ifndef DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
+#define DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
+
+#include <filesystem>
+
+
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/CorrectorComputer.hh>
+
+using namespace MatrixOperations;
+using std::shared_ptr;
+using std::make_shared;
+using std::string;
+using std::cout;
+using std::endl;
+
+// template <class Basis>
+// class EffectiveQuantitiesComputer : public CorrectorComputer<Basis> {
+
+template <class Basis>
+class EffectiveQuantitiesComputer {
+
+public:
+ static const int dimworld = 3;
+ // static const int nCells = 4;
+
+
+ static const int dim = Basis::GridView::dimension;
+
+ using Domain = typename CorrectorComputer<Basis>::Domain;
+
+ using VectorRT = typename CorrectorComputer<Basis>::VectorRT;
+ using MatrixRT = typename CorrectorComputer<Basis>::MatrixRT;
+
+ using Func2Tensor = typename CorrectorComputer<Basis>::Func2Tensor;
+ using FuncVector = typename CorrectorComputer<Basis>::FuncVector;
+
+ using VectorCT = typename CorrectorComputer<Basis>::VectorCT;
+
+ using HierarchicVectorView = typename CorrectorComputer<Basis>::HierarchicVectorView;
+
+protected:
+
+ CorrectorComputer<Basis>& correctorComputer_;
+ Func2Tensor prestrain_;
+
+public:
+ VectorCT B_load_TorusCV_; //<B, Chi>_L2
+ FieldMatrix<double, dim, dim> Q_; //effective moduli <LF_i, F_j>_L2
+ FieldVector<double, dim> Bhat_; //effective loads induced by prestrain <LF_i, B>_L2
+ FieldVector<double, dim> Beff_; //effective strains Mb = ak
+
+ // corrector parts
+ VectorCT phi_E_TorusCV_; //phi_i * (a,K)_i
+ VectorCT phi_perp_TorusCV_;
+ VectorCT phi_TorusCV_;
+ VectorCT phi_1_; //phi_i * (a,K)_i
+ VectorCT phi_2_;
+ VectorCT phi_3_;
+
+ // is this really interesting???
+ // double phi_E_L2norm_;
+ // double phi_E_H1seminorm_;
+
+ // double phi_perp_L2norm_;
+ // double phi_perp_H1seminorm_;
+
+ // double phi_L2norm_;
+ // double phi_H1seminorm_;
+
+ // double Chi_E_L2norm_;
+ // double Chi_perp_L2norm_;
+ // double Chi_L2norm_;
+
+
+ double B_energy_; // < B, B >_L B = F + Chi_perp + B_perp
+ double F_energy_; // < F, F >_L
+ double Chi_perp_energy_; // < Chi_perp, Chi_perp >_L
+ double B_perp_energy_; // < B_perp, B_perp >_L
+
+ //Chi(phi) is only implicit computed, can we store this?
+
+
+ ///////////////////////////////
+ // constructor
+ ///////////////////////////////
+ // EffectiveQuantitiesComputer(CorrectorComputer<Basis>& correctorComputer, Func2Tensor prestrain)
+ // : correctorComputer_(correctorComputer), prestrain_(prestrain)
+ EffectiveQuantitiesComputer(CorrectorComputer<Basis>& correctorComputer)
+ : correctorComputer_(correctorComputer)
+ {
+
+ // computePrestressLoadCV();
+ // computeEffectiveStrains();
+
+ // compute_phi_E_TorusCV();
+ // compute_phi_perp_TorusCV();
+ // compute_phi_TorusCV();
+
+ // computeCorrectorNorms();
+ // computeChiNorms();
+ // computeEnergiesPrestainParts();
+
+ // writeInLogfile();
+ }
+
+
+private:
+
+ void computeQ()
+ {
+
+ auto test = correctorComputer_.getLoad_alpha1();
+ auto phiContainer = correctorComputer_.getPhicontainer();
+
+ // std::cout << mu_ << "mu_ can be used this way"
+ // auto mu =
+ // auto lamba =
+ // auto gamma =
+
+ return ;
+ }
+
+
+
+
+
+
+
+
+
+ // void computePrestressLoadCV()
+ // {
+ // cout << "computePrestressLoadCV ..." << endl;
+ // auto cellAssembler = cellSolver_.getCellAssembler();
+ // cellAssembler.assembleLoadVector(prestrain_, B_load_TorusCV_);
+ // }
+
+ // void computeEffectiveStrains()
+ // {
+ // cout << "Compute effective strains ..." << endl;
+ // auto cellAssembler = cellSolver_.getCellAssembler();
+
+ // //E_h = (U + K e_cs, 0, 0)
+
+ // // lateral stretch strain E_U1 = e1 ot e1 = MatrixRT{{1, 0, 0}, {0, 0, 0}, {0, 0, 0}};
+ // Func2Tensor 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 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 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 E_K3 = [&] (const Domain& z) { return biotStrainApprox({0, 0, 0}, {0, 0, 1}, {0, z[dim-2], z[dim-1]}); };
+
+
+ // shared_ptr<VectorCT> phi_Ei_TorusCV[nCells] = {cellSolver_.getCorr_a(),
+ // cellSolver_.getCorr_K1(),
+ // cellSolver_.getCorr_K1(),
+ // cellSolver_.getCorr_K1()};
+
+ // shared_ptr<VectorCT> loadsCell[nCells] = {cellAssembler.getLoad_a(),
+ // cellAssembler.getLoad_K1(),
+ // cellAssembler.getLoad_K2(),
+ // cellAssembler.getLoad_K3()};
+
+ // const std::array<shared_ptr<Func2Tensor>, nCells> strainsE = {make_shared<Func2Tensor>(E_a),
+ // make_shared<Func2Tensor>(E_K1),
+ // make_shared<Func2Tensor>(E_K2),
+ // make_shared<Func2Tensor>(E_K3)};
+
+
+ // for (unsigned int k = 0; k < nCells; k++)
+ // for (unsigned int l = 0; l < nCells; l++){
+ // std::cout << "<L Fi,Fj> ... " << k << l << std::endl;
+ // M_[k][l] = cellAssembler.energyScalarProduct_FiFj(*strainsE[k], *phi_Ei_TorusCV[k], *strainsE[l], *phi_Ei_TorusCV[l]);
+ // }
+
+
+ // for (unsigned int k = 0; k < nCells; k++){
+ // std::cout << "<L Fi,B> ... " << k << std::endl;
+ // b_[k] = cellAssembler.energyScalarProduct_FiB(*strainsE[k], *phi_Ei_TorusCV[k], prestrain_);
+ // }
+
+
+ // M_.solve( aK_ ,b_);
+ // std::cout << "\nstretch twist bend1 bend2\naK = "<< aK_ << std::endl;
+ // }
+
+ // void compute_phi_E_TorusCV()
+ // {
+ // cout << "compute_phi_E_TorusCV ..." << endl;
+ // phi_E_TorusCV_ = *(cellSolver_.getCorr_a());
+ // phi_E_TorusCV_ *= aK_[0];
+ // auto tmp_xK1 = *(cellSolver_.getCorr_K1());
+ // auto tmp_xK2 = *(cellSolver_.getCorr_K2());
+ // auto tmp_xK3 = *(cellSolver_.getCorr_K3());
+ // tmp_xK1 *= aK_[1];
+ // tmp_xK2 *= aK_[2];
+ // tmp_xK3 *= aK_[3];
+ // phi_E_TorusCV_ += tmp_xK1;
+ // phi_E_TorusCV_ += tmp_xK2;
+ // phi_E_TorusCV_ += tmp_xK3;
+ // }
+
+ // void compute_phi_perp_TorusCV()
+ // {
+ // cout << "compute_phi_perp_TorusCV ..." << endl;
+ // cellSolver_.solvePrestrainCorrector(phi_perp_TorusCV_, B_load_TorusCV_);
+ // }
+
+ // void compute_phi_TorusCV()
+ // {
+ // cout << "compute_phi_TorusCV ..." << endl;
+ // phi_TorusCV_.resize(phi_E_TorusCV_.size());
+ // for (int i=0; i < phi_TorusCV_.size(); i++)
+ // phi_TorusCV_[i] = phi_E_TorusCV_[i] + phi_perp_TorusCV_[i];
+ // }
+
+
+
+
+
+
+
+
+ public:
+ ///////////////////////////////
+ // getter
+ ///////////////////////////////
+ CorrectorComputer<Basis> getCorrectorComputer(){return correctorComputer_;}
+
+ const shared_ptr<Basis> getBasis()
+ {
+ return correctorComputer_.getBasis();
+ }
+
+
+
+
+
+
+
+
+
+}; // end class
+
+
+
+
+
+#endif
\ No newline at end of file
diff --git a/src/Cell-Problem-New.cc b/src/Cell-Problem-New.cc
index 8c91ae68bd91d3aabef0d81892b6956b0c777498..e565d5f49da746bdea7a0250b629adf3d226e2a6 100644
--- a/src/Cell-Problem-New.cc
+++ b/src/Cell-Problem-New.cc
@@ -51,7 +51,8 @@
#include <dune/microstructure/prestrain_material_geometry.hh>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/vtk_filler.hh> //TEST
-#include <dune/microstructure/CorrectorComputer.hh> //TEST
+#include <dune/microstructure/CorrectorComputer.hh>
+#include <dune/microstructure/EffectiveQuantitiesComputer.hh>
#include <dune/solvers/solvers/umfpacksolver.hh> //TEST
@@ -1714,21 +1715,32 @@ int main(int argc, char *argv[])
// Compute reduced model
std::cout << "\ncompute effective model\n";
- auto cellAssembler = CellAssembler(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
+
+ // define type of FE-Basis...
+
- cellAssembler.solve();
- cellAssembler.computeNorms();
- cellAssembler.writeCorrectorsVTK(level);
- cellAssembler.check_Orthogonality();
+ // typedef Dune::Functions::LagrangeBasis<GridView, order> FEBasis;
+
+ // auto correctorComputer = CorrectorComputer(Basis_CE, muTerm, lambdaTerm, gamma, log, parameterSet);
+
+ // correctorComputer.solve();
+ // correctorComputer.computeNorms();
+ // correctorComputer.writeCorrectorsVTK(level);
+ // correctorComputer.check_Orthogonality();
+
+ // // -------------------------------------------
+ // auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer);
+ // // effectiveModelComputer.getQuantities(Q_eff,B_eff);
+ // effectiveQuantitiesComputer.computeQ();
+ std::cout << "\n WORKED \n";
// break;
- //TODO
- // cellAssembler.computeNorms();
- // cellAssembler.writeVTK();
+
+ // -------------------------------------------
- std::cout << "\n WORKED \n";
+
//TEST
// std::cout << "Test crossSectionDirectionScaling:" << std::endl;
@@ -2179,7 +2191,6 @@ int main(int argc, char *argv[])
std::cout << "---- (" << a << "," << b << ") ---- " << std::endl;
-
std::cout << "check_Orthogonality:" << check_Orthogonality(Basis_CE, muLocal, lambdaLocal,gamma,mContainer[a],mContainer[b],*phiDerContainer[a],*phiDerContainer[b],x3MatrixBasis[a]) << std::endl;