From 3df9afa1835cad8da1f6835f8e7290915c168cfa Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Thu, 19 Oct 2023 16:53:01 +0200
Subject: [PATCH] Cleanup
---
dune/microstructure/CorrectorComputer.hh | 421 +++---------------
.../EffectiveQuantitiesComputer.hh | 209 +--------
dune/microstructure/matrix_operations.hh | 179 +-------
dune/microstructure/prestrainedMaterial.hh | 340 ++------------
src/Cell-Problem.cc | 71 +--
5 files changed, 131 insertions(+), 1089 deletions(-)
diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 0c8dbcf6..d8591196 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -2,28 +2,27 @@
#define DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#include <map>
-
#include <dune/common/parametertree.hh>
-// #include <dune/common/float_cmp.hh>
-#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
-#include <dune/istl/matrixindexset.hh>
+
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
-#include <dune/microstructure/matrix_operations.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+
+#include <dune/istl/matrixindexset.hh>
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
-#include <dune/solvers/solvers/umfpacksolver.hh>
+#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/voigthelper.hh>
-using namespace MatrixOperations;
+#include <dune/solvers/solvers/umfpacksolver.hh>
+using namespace MatrixOperations;
using std::shared_ptr;
using std::make_shared;
using std::fstream;
-
template <class Basis, class Material> //, class LocalScalar, class Local2Tensor> // LocalFunction derived from basis?
class CorrectorComputer {
@@ -44,10 +43,8 @@ public:
using ElementMatrixCT = Dune::Matrix<double>;
protected:
-//private:
const Basis& basis_;
-
const Material& material_;
double gamma_;
@@ -55,9 +52,6 @@ protected:
fstream& log_; // Output-log
const Dune::ParameterTree& parameterSet_;
- // const FuncScalar mu_;
- // const FuncScalar lambda_;
-
MatrixCT stiffnessMatrix_;
VectorCT load_alpha1_,load_alpha2_,load_alpha3_; //right-hand side(load) vectors
@@ -75,7 +69,7 @@ protected:
const std::array<VoigtVector<double,3>,3 > matrixBasis_;
Func2Tensor x3G_1_ = [] (const Domain& x) {
- return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
+ return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
Func2Tensor x3G_2_ = [] (const Domain& x) {
@@ -109,12 +103,7 @@ public:
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))}),
phiOffset_(basis.size())
- {
-
- // assemble();
-
- //write VTK call here...
- }
+ {}
// -----------------------------------------------------------------
@@ -130,8 +119,6 @@ public:
assembleCellLoad(load_alpha3_ ,x3G_3_);
};
-
-
///////////////////////////////
// Getter
///////////////////////////////
@@ -147,26 +134,15 @@ public:
shared_ptr<VectorCT> getLoad_alpha1(){return make_shared<VectorCT>(load_alpha1_);}
shared_ptr<VectorCT> getLoad_alpha2(){return make_shared<VectorCT>(load_alpha2_);}
shared_ptr<VectorCT> getLoad_alpha3(){return make_shared<VectorCT>(load_alpha3_);}
-
- // shared_ptr<FuncScalar> getMu(){return make_shared<FuncScalar>(mu_);}
- // shared_ptr<FuncScalar> getLambda(){return make_shared<FuncScalar>(lambda_);}
-
shared_ptr<Material> getMaterial(){return make_shared<Material>(material_);}
// --- Get Correctors
- // shared_ptr<VectorCT> getMcontainer(){return make_shared<VectorCT>(mContainer);}
- // auto getMcontainer(){return make_shared<std::array<MatrixRT*, 3 > >(mContainer);}
auto getMcontainer(){return mContainer;}
shared_ptr<std::array<VectorCT, 3>> getPhicontainer(){return make_shared<std::array<VectorCT, 3>>(phiContainer);}
-
- // shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
auto getMatrixBasiscontainer(){return make_shared<std::array<VoigtVector<double,3>,3 >>(matrixBasis_);}
- // auto getx3MatrixBasiscontainer(){return make_shared<std::array<Func2Tensor, 3>>(x3MatrixBasisContainer_);}
auto getx3MatrixBasiscontainer(){return x3MatrixBasisContainer_;}
-
- // shared_ptr<VectorCT> getCorr_a(){return make_shared<VectorCT>(x_a_);}
shared_ptr<VectorCT> getCorr_phi1(){return make_shared<VectorCT>(phi_1_);}
shared_ptr<VectorCT> getCorr_phi2(){return make_shared<VectorCT>(phi_2_);}
shared_ptr<VectorCT> getCorr_phi3(){return make_shared<VectorCT>(phi_3_);}
@@ -174,12 +150,15 @@ public:
- // Get the occupation pattern of the stiffness matrix
+
+ /**
+ * @brief Get the occupation pattern of the stiffness matrix. Layout:
+ * | phi*phi m*phi |
+ * | phi *m m*m |
+ * @param nb MatrixIndexSet
+ */
void getOccupationPattern(Dune::MatrixIndexSet& nb)
{
- // OccupationPattern:
- // | phi*phi m*phi |
- // | phi *m m*m |
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
@@ -237,14 +216,11 @@ public:
std::cout << "finished occupation pattern\n";
}
- // template<class localFunction1, class localFunction2>
- // void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
- // ElementMatrixCT& elementMatrix,
- // const localFunction1& mu,
- // const localFunction2& lambda
- // )
- /** \brief Compute the stiffness matrix for one element
+ /** \brief Compute the stiffness matrix for one element. The local stiffness matrix
+ * has the form:
+ * | phi*phi m*phi |
+ * | phi *m m*m |
*
* \param quadRule The quadrature rule to be used
* \param phaseAtQuadPoint The material phase (i.e., the type of material) at each quadrature point
@@ -255,26 +231,18 @@ public:
ElementMatrixCT& elementMatrix
)
{
- // Local StiffnessMatrix of the form:
- // | phi*phi m*phi |
- // | phi *m m*m |
+
auto element = localView.element();
auto geometry = element.geometry();
- // using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
elementMatrix.setSize(localView.size()+3, localView.size()+3); //extend by dim ´R_sym^{2x2}
elementMatrix = 0;
-
- // auto elasticityTensor = material_.getElasticityTensor();
-
// LocalBasis-Offset
const int localPhiOffset = localView.size();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
- // std::cout << "localView.size(): " << localView.size() << std::endl;
- // std::cout << "nSf : " << nSf << std::endl;
int QPcounter= 0;
for (const auto& quadPoint : quadRule)
@@ -316,62 +284,7 @@ public:
for (size_t k=0; k < dimworld; k++)
for (size_t i=0; i < nSf; i++)
{
- // auto etmp = material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos));
- // auto etmp = elasticityTensor(defGradientU,element.geometry().global(quadPos));
- // auto etmp = material_.applyElasticityTensorLocal(defGradientU,quadPos);
- // printmatrix(std::cout, etmp, "etmp", "--");
- // double energyDensity= scalarProduct(etmp,defGradientV);
- // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),defGradientV);
-
- // auto tmp_1 = scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
- // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
- // if (abs(tmp_1-tmp_2)>1e-2)
- // {
- // std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
- // std::cout << "scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))<< std::endl;
- // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV)<< std::endl;
- // }
-
-
- // Python::Module module = Python::import("material_neukamm");
- // auto PythonindicatorFunction = Python::make_function<int>(module.get("f"));
- // if (PythonindicatorFunction(element.geometry().global(quadPos)) != material_.applyIndicatorFunction(element.geometry().global(quadPos)))
- // {
- // std::cout <<" element.geometry().global(quadPos): " << element.geometry().global(quadPos) << std::endl;
- // std::cout << "PythonindicatorFunction(element.geometry().global(quadPos)): " << PythonindicatorFunction(element.geometry().global(quadPos)) << std::endl;
- // // std::cout << "mu(quadPos):" << mu(quadPos) << std::endl;
- // std::cout << "indicatorFunction(element.geometry().global(quadPos)): " << indicatorFunction(element.geometry().global(quadPos)) << std::endl;
- // std::cout << "material_.applyIndicatorFunction(element.geometry().global(quadPos)): " << material_.applyIndicatorFunction(element.geometry().global(quadPos)) << std::endl;
- // std::cout << "indicatorFunction_material_1: " << indicatorFunction_material_1(element.geometry().global(quadPos)) << std::endl;
- // }
- // double energyDensity= scalarProduct(material_.ElasticityTensorGlobal(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
-
-
- // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
- // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
- // }
-
-
- // std::cout << "--------------- \n";
- // printmatrix(std::cout, defGradientU, "defGradientU", "--");
- // printmatrix(std::cout, material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.applyElasticityTensor(defGradientU,localIndicatorFunction(quadPos))", "--");
- // printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
-
-
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
- // double energyDensity= scalarProduct(material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(elasticityTensor(defGradientU,localIndicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
- // double energyDensity= scalarProduct(elasticityTensor(defGradientU,indicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,element.geometry().global(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.localElasticityTensor(defGradientU,quadPos,element),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.applyElasticityTensor(defGradientU,quadPos),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.applyElasticityTensorLocal(defGradientU,quadPos),defGradientV);
-
-
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV);
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), defGradientU, defGradientV); //TEST
- // double energyDensity = generalizedDensity(mu(quadPos), lambda(quadPos), defGradientU, defGradientV); // also works..
size_t col = localView.tree().child(k).localIndex(i);
@@ -381,92 +294,21 @@ public:
// "m*phi" & "phi*m" - part
for( size_t m=0; m<3; m++)
{
-
-
- // auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) ;
- // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
- // if (abs(tmp_1-tmp_2)>1e-2)
- // {
- // std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
- // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) " << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) << std::endl;
- // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
- // }
-
-
- // double energyDensityGphi = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
-
-
-
-
-
double energyDensityGphi = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),deformationGradient[j][l]);
- // double energyDensityGphi = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(defGradientV));
- // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV))" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV)) <<std::endl;
- // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV)" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV) << std::endl;
- // double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
- // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
- // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(defGradientV));
- // double energyDensityGphi= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(defGradientV));
- // double energyDensityGphi= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(defGradientV));
- // double energyDensityGphi= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),defGradientV);
- // double energyDensityGphi= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),defGradientV);
- // double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], defGradientV);
- // double energyDensityGphi = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], defGradientV); //TEST
+
auto value = energyDensityGphi * quadPoint.weight() * integrationElement;
elementMatrix[row][localPhiOffset+m] += value;
elementMatrix[localPhiOffset+m][row] += value;
}
-
}
// "m*m"-part
for(size_t m=0; m<3; m++) //TODO ist symmetric.. reicht die hälfte zu berechnen!!!
for(size_t n=0; n<3; n++)
{
-
- // std::cout << "QPcounter: " << QPcounter << std::endl;
- // std::cout << "m= " << m << " n = " << n << std::endl;
- // printmatrix(std::cout, basisContainer[m] , "basisContainer[m]", "--");
- // printmatrix(std::cout, basisContainer[n] , "basisContainer[n]", "--");
- // std::cout << "integrationElement:" << integrationElement << std::endl;
- // std::cout << "quadPoint.weight(): "<< quadPoint.weight() << std::endl;
- // std::cout << "mu(quadPos): " << mu(quadPos) << std::endl;
- // std::cout << "lambda(quadPos): " << lambda(quadPos) << std::endl;
-
- // auto tmp_1 = scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
- // auto tmp_2 = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
- // if (abs(tmp_1-tmp_2)>1e-2)
- // {
- // std::cout << "difference :" << abs(tmp_1-tmp_2) << std::endl;
- // std::cout << "scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])):" << scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n])) << std::endl;
- // std::cout << "linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])" << linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n])<< std::endl;
-
- // }
- // double energyDensityGG = scalarProduct(material_.ElasticityTensorGlobal(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
-
-
-
-
double energyDensityGG = voigtScalarProduct(material_.applyElasticityTensor(matrixBasis_[m],phase),matrixBasis_[n]);
- // double energyDensityGG = scalarProduct(material_.ElasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
- // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],localIndicatorFunction(quadPos)),sym(basisContainer[n]));
- // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[n]));
- // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
- // double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),basisContainer[n]);
- // double energyDensityGG= scalarProduct(elasticityTensor(basisContainer[m],indicatorFunction(quadPos)),sym(basisContainer[n]));
- // double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],element.geometry().global(quadPos)),sym(basisContainer[n]));
- // double energyDensityGG= scalarProduct(material_.applyElasticityTensor(basisContainer[m],quadPos),sym(basisContainer[n]));
- // double energyDensityGG= scalarProduct(material_.localElasticityTensor(basisContainer[m],quadPos,element),sym(basisContainer[n]));
- // double energyDensityGG= scalarProduct(material_.applyElasticityTensorLocal(basisContainer[m],quadPos),basisContainer[n]);
-
- // double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), basisContainer[m], basisContainer[n]);
- // double energyDensityGG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), basisContainer[m], basisContainer[n]); //TEST
+
elementMatrix[localPhiOffset+m][localPhiOffset+n] += energyDensityGG * quadPoint.weight() * integrationElement; // += !!!!! (Fixed-Bug)
- // std::cout << "energyDensityGG:" << energyDensityGG<< std::endl;
- // std::cout << "product " << energyDensityGG * quadPoint.weight() * integrationElement << std::endl;
- // printmatrix(std::cout, elementMatrix, "elementMatrix", "--");
}
}
// std::cout << "Number of QuadPoints:" << QPcounter << std::endl;
@@ -474,43 +316,28 @@ public:
}
-
- // 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
- // )
+ /**
+ * @brief Compute element load vector. Layout:
+ * | f*phi|
+ * | f*m |
+ *
+ * @tparam Vector
+ * @tparam LocalForce
+ * @param localView
+ * @param elementRhs
+ * @param forceTerm
+ */
template<class Vector, class LocalForce>
void computeElementLoadVector( const typename Basis::LocalView& localView,
Vector& elementRhs,
const LocalForce& forceTerm
)
{
- // | f*phi|
- // | --- |
- // | f*m |
- // using Element = typename LocalView::Element;
const auto element = localView.element();
const auto geometry = element.geometry();
- // constexpr int dim = Element::dimension;
- // constexpr int dimworld = dim;
- // using MatrixRT = FieldMatrix< double, dimworld, dimworld>;
-
- // auto elasticityTensor = material_.getElasticityTensor();
- // auto indicatorFunction = material_.getIndicatorFunction();
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
localIndicatorFunction.bind(element);
- // // auto indicatorFunction = material_.getLocalIndicatorFunction();
- // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
- // auto indicatorFunction = localFunction(indicatorFunctionGVF);
- // indicatorFunction.bind(element);
- // Func2int indicatorFunction = *material_.getIndicatorFunction();
-
// Set of shape functions for a single element
const auto& localFiniteElement= localView.tree().child(0).finiteElement();
@@ -522,14 +349,8 @@ public:
// LocalBasis-Offset
const int localPhiOffset = localView.size();
- // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
- // int orderQR = 0;
- // int orderQR = 1;
- // int orderQR = 2;
- // int orderQR = 3;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
- // std::cout << "Quad-Rule order used: " << orderQR << std::endl;
for (const auto& quadPoint : quad)
{
@@ -543,22 +364,6 @@ public:
for (size_t i=0; i< jacobians.size(); i++)
jacobians[i] = jacobians[i] * geometryJacobianInverse;
- //TEST
- // std::cout << "forceTerm(element.geometry().global(quadPos)):" << std::endl;
- // std::cout << forceTerm(element.geometry().global(quadPos)) << std::endl;
- // std::cout << "forceTerm(quadPos)" << std::endl;
- // std::cout << forceTerm(quadPos) << std::endl;
- //
- // //TEST QUadrature
- // std::cout << "quadPos:" << quadPos << std::endl;
- // std::cout << "element.geometry().global(quadPos):" << element.geometry().global(quadPos) << std::endl;
- // // //
- // // //
- // std::cout << "quadPoint.weight() :" << quadPoint.weight() << std::endl;
- // std::cout << "integrationElement(quadPos):" << integrationElement << std::endl;
- // std::cout << "mu(quadPos) :" << mu(quadPos) << std::endl;
- // std::cout << "lambda(quadPos) :" << lambda(quadPos) << std::endl;
-
// "f*phi"-part
for (size_t i=0; i < nSf; i++)
@@ -568,31 +373,12 @@ public:
MatrixRT tmpDefGradientV(0);
tmpDefGradientV[k][0] = jacobians[i][0][0]; // Y
tmpDefGradientV[k][1] = jacobians[i][0][1]; // X2
- // tmpDefGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2]; //
tmpDefGradientV[k][2] = jacobians[i][0][2]; // X3
VoigtVector<double,3> defGradientV = symVoigt(crossSectionDirectionScaling((1.0/gamma_),tmpDefGradientV));
- // double energyDensity= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
-
-
-
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),defGradientV);
- // double energyDensity= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(defGradientV));
- // double energyDensity= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(defGradientV));
- // double energyDensity= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),defGradientV);
- // double energyDensity= scalarProduct(material_.applyElasticityTensorLocal((-1.0)*forceTerm(quadPos),quadPos),defGradientV);
-
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV );
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)),forceTerm(quadPos), defGradientV ); //TEST
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),(-1.0)*forceTerm(quadPos), defGradientV ); //TEST
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos),forceTerm(element.geometry().global(quadPos)), defGradientV ); //TEST
-
+
size_t row = localView.tree().child(k).localIndex(i);
elementRhs[row] += energyDensity * quadPoint.weight() * integrationElement;
}
@@ -600,25 +386,9 @@ public:
// "f*m"-part
for (size_t m=0; m<3; m++)
{
- // double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
-
-
double energyDensityfG= voigtScalarProduct(material_.applyElasticityTensor((-1.0)*matrixToVoigt(forceTerm(quadPos)),localIndicatorFunction(quadPos)),matrixBasis_[m]);
- // double energyDensityfG= scalarProduct(material_.ElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
- // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(quadPos)),sym(basisContainer[m]));
- // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(element.geometry().global(quadPos))),sym(basisContainer[m]));
- // double energyDensityfG= scalarProduct(elasticityTensor((-1.0)*forceTerm(quadPos),indicatorFunction(quadPos)),sym(basisContainer[m]));
- // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
- // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),sym(basisContainer[m]));
- // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),basisContainer[m]);
- // double energyDensityfG = scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),quadPos),basisContainer[m]);
- // double energyDensityfG = scalarProduct(material_.localElasticityTensor((-1.0)*forceTerm(quadPos),quadPos,element),basisContainer[m]);
- // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] );
- // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(element.geometry().global(quadPos)), lambda(element.geometry().global(quadPos)), forceTerm(quadPos),basisContainer[m] ); //TEST
- // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), (-1.0)*forceTerm(quadPos),basisContainer[m] ); //TEST
- // double energyDensityfG = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), forceTerm(element.geometry().global(quadPos)),basisContainer[m] );//TEST
+
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
- // std::cout << "energyDensityfG * quadPoint.weight() * integrationElement: " << energyDensityfG * quadPoint.weight() * integrationElement << std::endl;
}
}
}
@@ -640,19 +410,11 @@ public:
bool cacheElementMatrices = true;
std::map<std::vector<int>, ElementMatrixCT> elementMatrixCache;
- // auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- // auto muLocal = localFunction(muGridF);
- // auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- // auto lambdaLocal = localFunction(lambdaGridF);
-
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
- // muLocal.bind(element);
- // lambdaLocal.bind(element);
+
const auto localPhiOffset = localView.size();
- // Dune::Timer Time;
- //std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
// Determine a suitable quadrature rule
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
@@ -697,14 +459,7 @@ public:
}
- // 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:
+ // TEST: Check Element-Stiffness-Symmetry:
if(parameterSet_.get<bool>("print_debug", false))
{
for (size_t i=0; i<localPhiOffset; i++)
@@ -755,31 +510,18 @@ public:
auto loadGVF = Dune::Functions::makeGridViewFunction(forceTerm, basis_.gridView());
auto loadFunctional = localFunction(loadGVF);
- // auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- // auto muLocal = localFunction(muGridF);
- // auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- // auto lambdaLocal = localFunction(lambdaGridF);
-
-
- // int counter = 1;
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
- // muLocal.bind(element);
- // lambdaLocal.bind(element);
loadFunctional.bind(element);
const auto localPhiOffset = localView.size();
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
VectorCT elementRhs;
- // std::cout << "----------------------------------Element : " << counter << std::endl; //TEST
- // counter++;
- // computeElementLoadVector(localView, muLocal, lambdaLocal, elementRhs, loadFunctional);
+
computeElementLoadVector(localView, elementRhs, loadFunctional);
- // computeElementLoadVector(localView, muLocal, lambdaLocal, gamma, elementRhs, forceTerm); //TEST
- // printvector(std::cout, elementRhs, "elementRhs", "--");
- // printvector(std::cout, elementRhs, "elementRhs", "--");
+
//////////////////////////////////////////////////////////////////////////////
// GLOBAL LOAD ASSEMBLY
//////////////////////////////////////////////////////////////////////////////
@@ -790,7 +532,6 @@ public:
}
for (size_t m=0; m<3; m++ )
b[phiOffset+m] += elementRhs[localPhiOffset+m];
- //printvector(std::cout, b, "b", "--");
}
}
@@ -818,8 +559,6 @@ public:
{
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++;
@@ -917,7 +656,6 @@ public:
r += localfun(quadPos) * quadPoint.weight() * integrationElement;
}
}
- // std::cout << "Domain-Area: " << area << std::endl;
return (1.0/area) * r ;
}
@@ -937,8 +675,17 @@ public:
}
- // -----------------------------------------------------------------
- // --- Solving the Corrector-problem:
+
+ /**
+ * @brief Solve corrector problem for different right-hand sides.
+ * Choose between solvers:
+ * 1 : CG-Solver
+ * 2 : GMRES
+ * 3 : QR (default)
+ * 4 : UMFPACK
+ *
+ * @return auto
+ */
auto solve()
{
bool set_oneBasisFunction_Zero = parameterSet_.get<bool>("set_oneBasisFunction_Zero", false);
@@ -961,13 +708,7 @@ public:
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);
@@ -1207,8 +948,11 @@ public:
}
- // -----------------------------------------------------------------
- // --- Write Correctos to VTK:
+ /**
+ * @brief Write correctors as grid functions to VTK.
+ *
+ * @param level GridLevel
+ */
void writeCorrectorsVTK(const int level)
{
std::string baseName = parameterSet_.get("baseName", "CellProblem-result");
@@ -1238,7 +982,6 @@ public:
{
computeL2Norm();
computeL2SymGrad();
-
std::cout<< "Frobenius-Norm of M1_: " << mContainer[0].frobenius_norm() << std::endl;
std::cout<< "Frobenius-Norm of M2_: " << mContainer[1].frobenius_norm() << std::endl;
std::cout<< "Frobenius-Norm of M3_: " << mContainer[2].frobenius_norm() << std::endl;
@@ -1334,10 +1077,12 @@ public:
return;
}
-
-
- // -----------------------------------------------------------------
- // --- Check Orthogonality relation Paper (75)
+ /**
+ * @brief Test: Check orthogonality relation (75) in
+ * [Böhnlein,Neukamm,Padilla-Garza,Sander - A homogenized bending theory for prestrained plates].
+ *
+ * @return auto
+ */
auto check_Orthogonality()
{
Dune::Timer orthogonalityTimer;
@@ -1355,13 +1100,6 @@ public:
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
-
-
- // auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
- // auto mu = localFunction(muGridF);
- // auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
- // auto lambda= localFunction(lambdaGridF);
-
for(int a=0; a<3; a++)
for(int b=0; b<3; b++)
{
@@ -1372,19 +1110,6 @@ public:
auto matrixFieldGGVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer_[a], basis_.gridView());
auto matrixFieldG = localFunction(matrixFieldGGVF);
- // auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
- // auto matrixFieldB = localFunction(matrixFieldBGVF);
-
- // using GridView = typename Basis::GridView;
- // using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- // using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
- // std::cout << "Press key.." << std::endl;
- // std::cin.get();
-
- // TEST
- // FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
- // printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
for (const auto& e : elements(basis_.gridView()))
{
@@ -1392,21 +1117,15 @@ public:
matrixFieldG.bind(e);
DerPhi1.bind(e);
DerPhi2.bind(e);
- // mu.bind(e);
- // lambda.bind(e);
localIndicatorFunction.bind(e);
double elementEnergy = 0.0;
- //double elementEnergy_HP = 0.0;
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
- // int orderQR = 0;
- // int orderQR = 1;
- // int orderQR = 2;
- // int orderQR = 3;
+
const auto& quad = Dune::QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -1414,27 +1133,20 @@ public:
const auto& quadPos = quadPoint.position();
const double integrationElement = geometry.integrationElement(quadPos);
-
auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + mContainer[b];
const auto strain1 = symVoigt(crossSectionDirectionScaling(1.0/gamma_, DerPhi1(quadPos)));
-
auto G = matrixFieldG(quadPos);
- // auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
-
+
auto tmp = matrixToVoigt(G + mContainer[a]) + strain1;
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),symVoigt(Chi));
- // double energyDensity= scalarProduct(material_.ElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), tmp, Chi);
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
- // elementEnergy += strain1 * quadPoint.weight() * integrationElement;
- //elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
+
}
energy += elementEnergy;
- //higherPrecEnergy += elementEnergy_HP;
}
if(parameterSet_.get<bool>("print_debug", false))
std::cout << "check_Orthogonality:" << "("<< a <<"," << b << "): " << energy << std::endl;
@@ -1487,7 +1199,6 @@ public:
}
-
}; // end class
#endif
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index 4deae7b2..b5fd0ef4 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -1,13 +1,12 @@
#ifndef DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_EFFECTIVEQUANTITIESCOMPUTER_HH
-// #include <filesystem>
-#include <dune/microstructure/matrix_operations.hh>
-#include <dune/microstructure/CorrectorComputer.hh>
+#include <dune/common/parametertree.hh>
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
#include <dune/istl/io.hh>
#include <dune/istl/matrix.hh>
-#include <dune/common/parametertree.hh>
+#include <dune/microstructure/matrix_operations.hh>
+#include <dune/microstructure/CorrectorComputer.hh>
using namespace MatrixOperations;
using std::shared_ptr;
@@ -16,8 +15,6 @@ using std::string;
using std::cout;
using std::endl;
-// template <class Basis>
-// class EffectiveQuantitiesComputer : public CorrectorComputer<Basis,Material> {
template <class Basis, class Material>
class EffectiveQuantitiesComputer {
@@ -36,10 +33,7 @@ public:
protected:
CorrectorComputer<Basis,Material>& correctorComputer_;
- // Func2Tensor prestrain_;
- // MatrixPhase prestrain_;
const Material& material_;
- // Func2int indicatorFunction_;
public:
@@ -55,10 +49,7 @@ public:
const Material& material)
: correctorComputer_(correctorComputer),
material_(material)
- {
-
- // prestrain_ = material_.getPrestrainFunction();
- }
+ {}
///////////////////////////////
@@ -77,14 +68,10 @@ public:
std::cout << "starting effective quantities computation..." << std::endl;
Dune::Timer effectiveQuantitiesTimer;
- // Get everything.. better TODO: with Inheritance?
- // auto test = correctorComputer_.getLoad_alpha1();
- // auto phiContainer = correctorComputer_.getPhicontainer();
auto MContainer = correctorComputer_.getMcontainer();
auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
- // auto mu_ = *correctorComputer_.getMu();
- // auto lambda_ = *correctorComputer_.getLambda();
+
auto gamma = correctorComputer_.getGamma();
auto basis = *correctorComputer_.getBasis();
Dune::ParameterTree parameterSet = correctorComputer_.getParameterSet();
@@ -94,12 +81,8 @@ public:
correctorComputer_.getCorr_phi3()
};
- // auto prestrainGVF = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
- // auto prestrainFunctional = localFunction(prestrainGVF);
-
auto localIndicatorFunction = material_.getLocalIndicatorFunction();
-
Qeff_ = 0 ;
Bhat_ = 0;
@@ -109,36 +92,22 @@ public:
double energy = 0.0;
double prestrain = 0.0;
auto localView = basis.localView();
- // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
+
auto GVFunc_a = derivative(Dune::Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
- // auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,phiContainer[b]));
auto localfun_a = localFunction(GVFunc_a);
- // auto localfun_b = localFunction(GVFunc_b);
-
- ///////////////////////////////////////////////////////////////////////////////
auto matrixFieldG1GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[a], basis.gridView());
auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
auto matrixFieldG2GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[b], basis.gridView());
auto matrixFieldG2 = localFunction(matrixFieldG2GVF);
- // auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
- // auto mu = localFunction(muGridF);
- // auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
- // auto lambda= localFunction(lambdaGridF);
-
- // using GridView = typename Basis::GridView;
-
for (const auto& element : elements(basis.gridView()))
{
localView.bind(element);
matrixFieldG1.bind(element);
matrixFieldG2.bind(element);
localfun_a.bind(element);
- // DerPhi2.bind(e);
- // mu.bind(e);
- // lambda.bind(e);
- // prestrainFunctional.bind(e);
+
localIndicatorFunction.bind(element);
double elementEnergy = 0.0;
@@ -147,12 +116,8 @@ public:
auto geometry = element.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
- // int orderQR = 0;
- // int orderQR = 1;
- // int orderQR = 2;
- // int orderQR = 3;
+
const auto& quad = Dune::QuadratureRules<double, dim>::rule(element.type(), orderQR);
for (const auto& quadPoint : quad)
@@ -162,35 +127,20 @@ public:
auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + MContainer[a];
-
auto G1 = matrixFieldG1(quadPos);
auto G2 = matrixFieldG2(quadPos);
- // auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
- // auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
auto X1 = matrixToVoigt(G1 + Chi1);
- // auto X2 = G2 + Chi2;
double energyDensity= voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)), matrixToVoigt(sym(G2)));
- // double energyDensity= scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, G2);
+
elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
if (b==0)
{
- // std::cout << "WORKS TILL HERE" << std::endl;
- // auto prestrainPhaseValue_old = prestrain_(localIndicatorFunction(quadPos));
auto prestrainPhaseValue = material_.prestrain(localIndicatorFunction(quadPos),element.geometry().global(quadPos));
- // if (localIndicatorFunction(quadPos) != 3)
- // {
- // std::cout << "element.geometry().global(quadPos):"<< element.geometry().global(quadPos) << std::endl;
- // printmatrix(std::cout, prestrainPhaseValue_old, "prestrainPhaseValue_old", "--");
- // printmatrix(std::cout, prestrainPhaseValue, "prestrainPhaseValue", "--");
- // }
auto value = voigtScalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),matrixToVoigt(sym(prestrainPhaseValue)));
elementPrestrain += value * quadPoint.weight() * integrationElement;
- // elementPrestrain += scalarProduct(material_.ElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(prestrainPhaseValue)) * quadPoint.weight() * integrationElement;
- // elementPrestrain += linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), X1, prestrainFunctional(quadPos)) * quadPoint.weight() * integrationElement;
}
}
energy += elementEnergy;
@@ -211,7 +161,7 @@ public:
// Compute effective Prestrain B_eff (by solving linear system)
//////////////////////////////
- // std::cout << "------- Information about Q matrix -----" << std::endl; // TODO
+ // std::cout << "------- Information about Q matrix -----" << std::endl; // TODO
// MatrixInfo<MatrixRT> matrixInfo(Qeff_,true,2,1);
// std::cout << "----------------------------------------" << std::endl;
Qeff_.solve(Beff_,Bhat_);
@@ -228,7 +178,6 @@ public:
log << "Beff_: " << Beff_ << " (Effective Prestrain)" << std::endl;
log << "------------------------ " << std::endl;
- // TEST
// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
std::cout << "Time required to solve for effective quantities: " << effectiveQuantitiesTimer.elapsed() << std::endl;
return ;
@@ -236,11 +185,10 @@ public:
// -----------------------------------------------------------------
- // --- write Data to Matlab / Optimization-Code
+ // --- write Data to .Txt-File for Matlab / Optimization-Code
void writeToMatlab(std::string outputPath)
{
std::cout << "write effective quantities as .txt files to output folder..." << std::endl;
- //writeMatrixToMatlab(Q, "../../Matlab-Programs/QMatrix.txt");
writeMatrixToMatlab(Qeff_, outputPath + "/QMatrix.txt");
// write effective Prestrain in Matrix for Output
Dune::FieldMatrix<double,1,3> BeffMat;
@@ -266,17 +214,12 @@ public:
auto matrixFieldA = localFunction(matrixFieldAGVF);
auto matrixFieldBGVF = Dune::Functions::makeGridViewFunction(matrixFieldFuncB, basis.gridView());
auto matrixFieldB = localFunction(matrixFieldBGVF);
- // auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
- // auto mu = localFunction(muGridF);
- // auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
- // auto lambda= localFunction(lambdaGridF);
+
for (const auto& e : elements(basis.gridView()))
{
localView.bind(e);
matrixFieldA.bind(e);
matrixFieldB.bind(e);
- // mu.bind(e);
- // lambda.bind(e);
localIndicatorFunction.bind(e);
double elementEnergy = 0.0;
@@ -292,8 +235,7 @@ public:
const double integrationElement = geometry.integrationElement(quadPos);
double energyDensity= scalarProduct(material_.applyElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
- // double energyDensity= scalarProduct(material_.ElasticityTensor(matrixFieldA(quadPos),localIndicatorFunction(quadPos)),sym(matrixFieldB(quadPos)));
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), matrixFieldA(quadPos), matrixFieldB(quadPos));
+
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
@@ -301,131 +243,6 @@ public:
return energy;
}
-
-
- // --- Alternative that does not use orthogonality relation (75) in the paper
- // void computeFullQ()
- // {
- // auto MContainer = correctorComputer_.getMcontainer();
- // auto MatrixBasisContainer = correctorComputer_.getMatrixBasiscontainer();
- // auto x3MatrixBasisContainer = correctorComputer_.getx3MatrixBasiscontainer();
- // auto mu_ = *correctorComputer_.getMu();
- // auto lambda_ = *correctorComputer_.getLambda();
- // auto gamma = correctorComputer_.getGamma();
- // auto basis = *correctorComputer_.getBasis();
-
- // shared_ptr<VectorCT> phiBasis[3] = {correctorComputer_.getCorr_phi1(),
- // correctorComputer_.getCorr_phi2(),
- // correctorComputer_.getCorr_phi3()
- // };
-
- // auto prestrainGVF = Dune::Functions::makeGridViewFunction(prestrain_, basis.gridView());
- // auto prestrainFunctional = localFunction(prestrainGVF);
-
- // Qeff_ = 0 ;
- // Bhat_ = 0;
-
- // for(size_t a=0; a < 3; a++)
- // for(size_t b=0; b < 3; b++)
- // {
- // double energy = 0.0;
- // double prestrain = 0.0;
- // auto localView = basis.localView();
- // // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiContainer[a]));
- // auto GVFunc_a = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[a]));
- // auto GVFunc_b = derivative(Functions::makeDiscreteGlobalBasisFunction<VectorRT>(basis,*phiBasis[b]));
- // auto localfun_a = localFunction(GVFunc_a);
- // auto localfun_b = localFunction(GVFunc_b);
-
- // ///////////////////////////////////////////////////////////////////////////////
- // auto matrixFieldG1GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[a], basis.gridView());
- // auto matrixFieldG1 = localFunction(matrixFieldG1GVF);
- // auto matrixFieldG2GVF = Dune::Functions::makeGridViewFunction(x3MatrixBasisContainer[b], basis.gridView());
- // auto matrixFieldG2 = localFunction(matrixFieldG2GVF);
-
- // auto muGridF = Dune::Functions::makeGridViewFunction(mu_, basis.gridView());
- // auto mu = localFunction(muGridF);
- // auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambda_, basis.gridView());
- // auto lambda= localFunction(lambdaGridF);
-
- // // using GridView = typename Basis::GridView;
-
- // for (const auto& e : elements(basis.gridView()))
- // {
- // localView.bind(e);
- // matrixFieldG1.bind(e);
- // matrixFieldG2.bind(e);
- // localfun_a.bind(e);
- // localfun_b.bind(e);
- // mu.bind(e);
- // lambda.bind(e);
- // prestrainFunctional.bind(e);
-
- // double elementEnergy = 0.0;
- // double elementPrestrain = 0.0;
-
- // auto geometry = e.geometry();
- // const auto& localFiniteElement = localView.tree().child(0).finiteElement();
-
- // // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 + 5 ); // TEST
- // int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
- // // 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 Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a] + matrixFieldG1(quadPos);
- // auto Chi2 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_b(quadPos))) + *MContainer[b] + matrixFieldG2(quadPos);
-
- // // auto G1 = matrixFieldG1(quadPos);
- // // auto G2 = matrixFieldG2(quadPos);
- // // auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
- // // auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
-
- // // auto X1 = G1 + Chi1;
- // // auto X2 = G2 + Chi2;
-
-
- // double energyDensity = linearizedStVenantKirchhoffDensity(mu(quadPos), lambda(quadPos), Chi1, Chi2);
- // elementEnergy += energyDensity * quadPoint.weight() * integrationElement; // quad[quadPoint].weight() ???
- // }
- // energy += elementEnergy;
- // prestrain += elementPrestrain;
-
- // }
- // Qeff_[a][b] = energy;
- // if (b==0)
- // Bhat_[a] = prestrain;
- // }
- // printmatrix(std::cout, Qeff_, "Matrix Qeff_", "--");
- // printvector(std::cout, Bhat_, "Bhat_", "--");
- // ///////////////////////////////
- // // Compute effective Prestrain B_eff (by solving linear system)
- // //////////////////////////////
- // // std::cout << "------- Information about Q matrix -----" << std::endl; // TODO
- // // MatrixInfo<MatrixRT> matrixInfo(Qeff_,true,2,1);
- // // std::cout << "----------------------------------------" << std::endl;
- // Qeff_.solve(Beff_,Bhat_);
- // printvector(std::cout, Beff_, "Beff_", "--");
-
- // //LOG-Output
- // auto& log = *(correctorComputer_.getLog());
- // log << "--- Prestrain Output --- " << std::endl;
- // log << "Bhat_: " << Bhat_ << std::endl;
- // log << "Beff_: " << Beff_ << " (Effective Prestrain)" << std::endl;
- // log << "------------------------ " << std::endl;
- // return ;
- // }
-
-
}; // end class
-
-
#endif
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index ab0a5933..4e0a89c0 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -9,14 +9,6 @@ namespace MatrixOperations {
using std::sin;
using std::cos;
- // static MatrixRT sym (MatrixRT M) { // 1/2 (M^T + M)
- // MatrixRT ret(0);
- // for (int i = 0; i< 3; i++)
- // for (int j = 0; j< 3; j++)
- // ret[i][j] = 0.5 * (M[i][j] + M[j][i]);
- // return ret;
- // }
-
static MatrixRT sym (MatrixRT M) { // 1/2 (M^T + M)
MatrixRT ret(0);
@@ -62,9 +54,6 @@ namespace MatrixOperations {
return rankoneTensorproduct(U_plus_k_x_ecs, e_1);
}
-// static MatrixRT crossSectionDirectionScaling(double w, MatrixRT M){
-// return {M[0], M[1], w*M[2]};
-// }
static MatrixRT crossSectionDirectionScaling(double w, MatrixRT M){
return {{M[0][0], M[0][1], w*M[0][2]},
@@ -159,50 +148,17 @@ namespace MatrixOperations {
}
}
-// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){ // ?? Check with Robert
-// E1= sym(E1);
-// E2 = sym(E2);
-// double t1 = scalarProduct(E1,E2);
-// t1 *= 2* mu;
-// double t2 = trace(E1)*trace(E2);
-// t2 *= lambda;
-// return t1 + t2;
-// }
-//
-// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
-// {
-// auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
-// return scalarProduct(t1,sym(E2));
-// }
-// //
static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
{
auto t1 = 2.0 * mu * sym(E1) + MatrixRT(Dune::ScaledIdentityMatrix<double,3>(lambda * trace(sym(E1))));
auto tmp1 = scalarProduct(t1,sym(E2));
- // auto tmp1 = scalarProduct(t1,E2);
-// auto t2 = 2.0 * mu * sym(E2) + lambda * trace(sym(E2)) * Id();
-// auto tmp2 = scalarProduct(t2,sym(E1));
-//
-// if(abs(tmp1-tmp2) > 1e-8 )
-// {
-// std::cout << "linearizedStVenantKirchhoffDensity NOT Symmetric!" << std::endl;
-// std::cout << "tmp1: " << tmp1 << std::endl;
-// std::cout << "tmp2: " << tmp2 << std::endl;
-// }
return tmp1;
}
- // static MatrixRT elasticityTensor()
- // {
-
-
- // }
-
// --- Generalization: Define Quadratic QuadraticForm
-
static double QuadraticForm(const double mu, const double lambda, const MatrixRT M){
auto tmp1 = sym(M);
@@ -213,14 +169,7 @@ namespace MatrixOperations {
}
-// static double linearizedStVenantKirchhoffDensity(const double mu, const double lambda, MatrixRT F, MatrixRT G){
-// /// Write this whole File as a Class that uses lambda,mu as members ?
-//
-// // Define L via Polarization-Identity from QuadratifForm
-// // <LF,G> := (1/2)*(Q(F+G) - Q(F) - Q(G) )
-// return (1.0/2.0)*(QuadraticForm(mu,lambda,F+G) - QuadraticForm(mu,lambda,F) - QuadraticForm(mu,lambda,G) );
-// }
-
+
static double generalizedDensity(const double mu, const double lambda, MatrixRT F, MatrixRT G){
/// Write this whole File as a Class that uses lambda,mu as members ?
@@ -247,8 +196,6 @@ namespace MatrixOperations {
}
-
-
extern "C"
{
@@ -257,133 +204,9 @@ namespace MatrixOperations {
return sym(M);
}
-
-
-
-
}
-
-
-
-
-
- /*
- template<double phi>
- static bool isInRotatedPlane(double x1, double x2){
- return cos(phi)*x1 + sin(phi)*x2 > 0;
- }*/
-
}
-
-
-// OPTIONAL H1-Seminorm ...
-
-/*
-template<class VectorCT>
-double semi_H1_vectorScalarProduct(VectorCT& coeffVector1, VectorCT& coeffVector2)
-{
-
-double ret(0);
-auto localView = basis_.localView();
-
-for (const auto& e : elements(basis_.gridView()))
- {
- localView.bind(e);
-
- auto geometry = e.geometry();
-
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- const auto nSf = localFiniteElement.localBasis().size();
-
- int order = 2*(dim*localFiniteElement.localBasis().order()-1);
- const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), order);
-
- for (const auto& quadPoint : quad)
- {
- double elementScp(0);
-
- auto pos = quadPoint.position();
- const auto jacobian = geometry.jacobianInverseTransposed(pos);
- const double integrationElement = geometry.integrationElement(pos);
-
- std::vector<FieldMatrix<double,1,dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(pos, referenceGradients);
-
- std::vector< VectorRT > gradients(referenceGradients.size());
- for (size_t i=0; i< gradients.size(); i++)
- jacobian.mv(referenceGradients[i][0], gradients[i]);
-
- for (size_t i=0; i < nSf; i++)
- for (size_t j=0; j < nSf; j++)
- for (size_t k=0; k < dimworld; k++)
- for (size_t l=0; l < dimworld; l++)
- {
- MatrixRT defGradient1(0);
- defGradient1[k] = gradients[i];
-
- MatrixRT defGradient2(0);
- defGradient2[l] = gradients[j];
-
- size_t localIdx1 = localView.tree().child(k).localIndex(i);
- size_t globalIdx1 = localView.index(localIdx1);
-
- size_t localIdx2 = localView.tree().child(l).localIndex(j);
- size_t globalIdx2 = localView.index(localIdx2);
-
- double scp(0);
- for (int ii=0; ii<dimworld; ii++)
- for (int jj=0; jj<dimworld; jj++)
- scp += coeffVector1[globalIdx1] * defGradient1[ii][jj] * coeffVector2[globalIdx2] * defGradient2[ii][jj];
-
- elementScp += scp;
- }
-
- ret += elementScp * quadPoint.weight() * integrationElement;
- } //qp
- } //e
- return ret;
-}*/
-
-
-
-
-
-
-/*
- class BiotStrainApprox
- {
-
- // E_h = h^{-1} (R^T F_h - Id)
- //
- //E_h = (U + K e_cs, 0, 0)
- //
-
- Func2Tensor E_a = [] (const Domain& z) //extra Klasse
- {
- return biotStrainApprox({1, 0, 0}, {0, 0, 0}, {0, z[1], z[2]}); //MatrixRT{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
-
- // twist in x2-x3 plane E_K1 = sym (e1 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K1 = [] (const Domain& z)
- {
- return biotStrainApprox({0, 0, 0}, {1, 0, 0}, {0, z[1], z[2]}); //MatrixRT{{0.0,-0.5*z[2], 0.5*z[1]}, {-0.5*z[2], 0.0 , 0.0 }, {0.5*z[1],0.0,0.0}};
- };
-
- // bend strain in x1-x2 plane E_K2 = sym (e2 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K2 = [] (const Domain& z)
- {
- return biotStrainApprox({0, 0, 0}, {0, 1, 0}, {0, z[1], z[2]}); //MatrixRT{{z[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0}};
- };
-
- // bend strain in x1-x3 plane E_K3 = sym (e3 x (0,x2,x3)^T ot e1)
- Func2Tensor E_K3 = [] (const Domain& z)
- {
- return biotStrainApprox({0, 0, 0}, {0, 0, 1}, {0, z[1], z[2]}); //MatrixRT{{-z[1], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
- };
- };
-*/
-
#endif
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 2cc7fcfe..85af8567 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -1,19 +1,16 @@
#ifndef DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
#define DUNE_MICROSTRUCTURE_PRESTRAINEDMATERIAL_HH
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
-#include <dune/grid/yaspgrid.hh>
+#include <dune/common/parametertree.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
#include <dune/functions/gridfunctions/gridviewentityset.hh>
-#include <dune/common/parametertree.hh>
#include <dune/fufem/dunepython.hh>
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
+#include <dune/grid/yaspgrid.hh>
#include <dune/microstructure/matrix_operations.hh>
#include <dune/microstructure/voigthelper.hh>
-// #include <dune/microstructure/materialDefinitions.hh> //deprecated
-
-
using namespace MatrixOperations;
using namespace Dune::Functions;
using namespace Dune::Functions::BasisFactory;
@@ -22,61 +19,9 @@ using std::abs;
using std::sqrt;
using std::sin;
using std::cos;
-
using std::shared_ptr;
using std::make_shared;
-
-
-
-
-
-// // ----------------------------------------------------------------------------------
-// template<class Domain>
-// int indicatorFunction_material(const Domain& x)
-// {
-// double theta=0.25;
-// if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
-// return 1; //#Phase1
-// else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
-// return 2; //#Phase2
-// else
-// return 3; //#Phase3
-// }
-// MatrixRT getL(const int& phase)
-// {
-// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
-// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
-
-
-
-// if (phase == 1)
-// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id(); //#Phase1 //Isotrop(mu,lambda)
-// else if (phase == 2)
-// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id(); //#Phase2
-// else
-// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id(); //#Phase3
-// }
-// MatrixRT prestrain_material(const int& phase)
-// {
-// if (phase == 1)
-// return {{1.0, 0.0, 0.0}, //#Phase1
-// {0.0, 1.0, 0.0},
-// {0.0, 0.0, 1.0}
-// };
-// else if (phase == 2)
-// return {{1.0, 0.0, 0.0}, //#Phase2
-// {0.0, 1.0, 0.0},
-// {0.0, 0.0, 1.0}
-// };
-// else
-// return {{0.0, 0.0, 0.0}, //#Phase3
-// {0.0, 0.0, 0.0},
-// {0.0, 0.0, 0.0}
-// };
-// }
-
-
//--------------------------------------------------------
template <class GridView>
class prestrainedMaterial
@@ -102,7 +47,6 @@ public:
protected:
const GridView& gridView_;
const Dune::ParameterTree parameterSet_;
- // std::string materialFunctionName_;
// Elasticity tensors (in Voigt notation)
std::vector<VoigtMatrix<double,3> > L_;
@@ -119,33 +63,7 @@ protected:
LocalFunction<int(const Domain&), typename GridViewEntitySet<GridView, 0>::Element > localIndicatorFunction_;
- //Transformation matrix for Voigt notation
- Dune::FieldMatrix<double,6,6> D_ = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
- {0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
- {0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
- {0.0, 0.0, 0.0, sqrt(2.0) , 0.0, 0.0},
- {0.0, 0.0, 0.0, 0.0 , sqrt(2.0), 0.0},
- {0.0, 0.0, 0.0, 0.0 , 0.0, sqrt(2.0)}
- };
- Dune::FieldMatrix<double,6,6> D_inv = {{1.0, 0.0, 0.0, 0.0 , 0.0, 0.0},
- {0.0, 1.0, 0.0, 0.0 , 0.0, 0.0},
- {0.0, 0.0, 1.0, 0.0 , 0.0, 0.0},
- {0.0, 0.0, 0.0, 1.0/sqrt(2.0) , 0.0, 0.0},
- {0.0, 0.0, 0.0, 0.0 , 1.0/sqrt(2.0), 0.0},
- {0.0, 0.0, 0.0, 0.0 , 0.0, 1.0/sqrt(2.0)}
- };
-
-
-
public:
- ///////////////////////////////
- // constructor
- ///////////////////////////////
- // prestrainedMaterial(const GridView gridView,
- // const Dune::ParameterTree& parameterSet)
- // : gridView_(gridView),
- // parameterSet_(parameterSet)
- // {
prestrainedMaterial(const GridView gridView,
const Dune::ParameterTree& parameterSet,
Python::Module pyModule)
@@ -200,7 +118,7 @@ Func2Tensor setupPrestrainPhase(const int phase){
Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
{
std::string phaseType;
- // module_.get("phase" + std::to_string(phase) + "_type").toC<std::string>(phaseType);
+
phaseType = parameterSet_.get<std::string>("phase" + std::to_string(phase) + "_type");
std::cout << "Setup phase "<< phase << " as " << phaseType << " material." << std::endl;
@@ -213,15 +131,7 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
return isotropic(materialParameters[0],materialParameters[1]);
}
- /**
- * @brief For other material classes. The stiffness-/compliance matrix depends on the
- * coordinate frame.
- *
- * The default frame vectors:
- */
- // VectorRT frameVector1 = {1.0,0.0,0.0};
- // VectorRT frameVector2 = {0.0,1.0,0.0};
- // VectorRT frameVector3 = {0.0,0.0,1.0};
+
/**
* @brief For other material classes. The stiffness-/compliance matrix depends on the
* coordinate frame.
@@ -233,20 +143,11 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
try
{
- // module.get("phase" + std::to_string(phase) + "_FrameVector1").toC<VectorRT>(frameVector1);
- // module.get("phase" + std::to_string(phase) + "_FrameVector2").toC<VectorRT>(frameVector2);
- // module.get("phase" + std::to_string(phase) + "_FrameVector3").toC<VectorRT>(frameVector3);
- // printvector(std::cout, frameVector1, "frameVector1: ", "--");
- // printvector(std::cout, frameVector2, "frameVector2: ", "--");
- // printvector(std::cout, frameVector3, "frameVector3: ", "--");
- // module_.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
- // module_.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
axis = parameterSet_.get<int>("phase" + std::to_string(phase) + "_axis", 0);
angle = parameterSet_.get<double>("phase" + std::to_string(phase) + "_angle", 0);
}
catch(Dune::Exception& e)
{
- // std::cout << "(Could not read frame vectors for material phase) " << phase << ": - Canonical Frame (default) is used." << std::endl;
std::cout << "(Could not read rotation axis & angle for material phase) " << phase << ": - Canonical Frame (default) is used." << std::endl;
}
@@ -255,29 +156,24 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
{
Dune::FieldVector<double,9> materialParameters(0);
module_.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,9>>(materialParameters);
- // return orthotropic(frameVector1,frameVector2,frameVector3,materialParameters);
return orthotropic(axis,angle,materialParameters);
}
if(phaseType == "transversely_isotropic")
{
Dune::FieldVector<double,5> materialParameters(0);
module_.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldVector<double,5>>(materialParameters);
- // return transversely_isotropic(frameVector1,frameVector2,frameVector3,materialParameters);
return transversely_isotropic(axis,angle,materialParameters);
}
if(phaseType == "general_anisotropic")
{
Dune::FieldMatrix<double,6,6> materialParameters(0); //compliance matric
module_.get("materialParameters_phase" + std::to_string(phase)).toC<Dune::FieldMatrix<double,6,6>>(materialParameters);
- // return general_anisotropic(frameVector1,frameVector2,frameVector3,materialParameters);
return general_anisotropic(axis,angle,materialParameters);
}
else
DUNE_THROW(Dune::Exception, " Unknown material class for phase " + std::to_string(phase));
}
-
-
/**
* @brief Setup method that determines the correct stiffness-matrix (constitutive law)
* for each material phase separately.
@@ -290,7 +186,6 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
int phases;
- // module_.get("Phases").toC<int>(phases);
phases = parameterSet_.get<int>("Phases");
std::cout << "---- Starting material setup... (Number of Phases: "<< phases << ") " << std::endl;
@@ -306,28 +201,6 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
std::cout << "Material setup done." << std::endl;
}
- // void setup(const std::string name) // can't use materialFunctionName_ here?
- // {
- // Python::Module module = Python::import(name);
- // indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
-
- // int phases;
- // module.get("Phases").toC<int>(phases);
- // std::cout << "---- Starting material setup... (Number of Phases: "<< phases << ") " << std::endl;
-
- // L_.resize(phases);
- // prestrain_.resize(phases);
-
- // for (int i=0; i<phases; i++)
- // {
- // L_[i] = setupPhase(module,i+1);
- // prestrain_[i] = setupPrestrainPhase(module, i+1);
- // }
-
- // std::cout << "Material setup done." << std::endl;
- // }
-
-
////////////////////////////////////////////////////////
// MATERIAL CLASSES DEFINITIONS:
////////////////////////////////////////////////////////
@@ -345,14 +218,6 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
}
//-------------------------------------------------------------------------------
- //--- Definition of orthotropic elasticity tensor (in voigt notation)
- // - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
- // - Output: compliance Matrix
- // Dune::FieldMatrix<double,6,6> orthotropic(const VectorRT& framevector1,
- // const VectorRT& framevector2,
- // const VectorRT& framevector3,
- // const Dune::FieldVector<double,9>& p //elastic moduli {E1,E2,E3,G12,G23,G31,nu12,nu13,nu23}
- // )
/**
* @brief Definition of orthotropic elasticity/stiffness tensor (in Voigt notation).
@@ -407,7 +272,6 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
* @brief (optional) Tramsform compliance matrix to the correct frame.
*/
if(abs(angle) > 1e-12){
- // std::cout << "material frame rotated around axis: " << axis << " by an angle of: " << angle << std::endl;
std::cout << "material frame rotated around axis: " << axis_string << " by an angle of: " << angle << std::endl;
Dune::FieldMatrix<double,6,6> C_rot = rotationMatrixCompliance(axis,angle);
(Compliance.leftmultiply(C_rot)).rightmultiply(C_rot.transposed());
@@ -425,15 +289,14 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
}
//-------------------------------------------------------------------------------
-
- // Dune::FieldMatrix<double,6,6> transversely_isotropic(const VectorRT& framevector1,
- // const VectorRT& framevector2,
- // const VectorRT& framevector3,
- // const Dune::FieldVector<double,5>& p //elastic moduli {E1,E2,G12,nu12,nu23}
- // )
- //--- Definition of transversely isotropic elasticity tensor (in voigt notation)
- // - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
- // - Output: compliance Matrix
+ /**
+ * @brief Definition of transversely isotropic elasticity tensor (in Voigt notation)
+ * (see https://en.wikipedia.org/wiki/Poisson%27s_ratio)
+ * @param axis rotate material frame around axis
+ * @param angle rotation angle
+ * @param p elastic constants: {E1,E2,G12,nu12,nu23}
+ * @return Stiffness matrix
+ */
Dune::FieldMatrix<double,6,6> transversely_isotropic(const int& axis,
const double& angle,
const Dune::FieldVector<double,5>& p
@@ -483,9 +346,8 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
}
printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
- // printmatrix(std::cout, C, "C", "--");
- //--- return elasticity tensor
+ // invert Compliance matrix to obtain elasticity/stiffness matrix
Compliance.invert();
// printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
@@ -496,20 +358,19 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
//-------------------------------------------------------------------------------
- // --- Definition of transversely isotropic elasticity tensor (in voigt notation)
- // - Input: (Youngs modulus, shear modulus, Poisson ratio): {E1,E2,G12,nu12,nu23}
- // - Output: inverse compliance Matrix
- // Dune::FieldMatrix<double,6,6> general_anisotropic(const VectorRT& framevector1,
- // const VectorRT& framevector2,
- // const VectorRT& framevector3,
- // const Dune::FieldMatrix<double,6,6>& C //compliance matrix
- // )
+ /**
+ * @brief Definition of general anisotropic elasticity tensor (in Voigt notation)
+ *
+ * @param axis rotate material frame around axis
+ * @param angle rotation angle
+ * @param Compliance Compliance matrix
+ * @return Stiffness matrix
+ */
Dune::FieldMatrix<double,6,6> general_anisotropic(const int& axis,
const double& angle,
Dune::FieldMatrix<double,6,6>& Compliance
)
{
-
std::string axis_string;
switch (axis)
{
@@ -531,9 +392,7 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
-
-
- //--- return stiffness/elasticity tensor (in Voigt notation)
+ // invert Compliance matrix to obtain elasticity/stiffness matrix
Compliance.invert();
// printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
@@ -548,25 +407,24 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
{
VoigtVector<double,3> out(0);
- // printmatrix(std::cout, L_[0], "L_[0]", "--");
L_[phase-1].mv(G,out);
return out;
}
- ////////////////////////////////////////////////////////////////////////////////////////
MatrixRT prestrain(const int& phase,const Domain& x) const
{
return prestrain_[phase-1](x);
}
-
-////////////////////////////////////////////////////////
-
-
- //--- Apply elasticityTensor_ to input Matrix G at material phase
- // input: Matrix G, material-phase
+ /**
+ * @brief Apply elasticityTensor_ to input Matrix G at material phase
+ *
+ * @param G input matrix
+ * @param phase material phase
+ * @return MatrixRT
+ */
MatrixRT ElasticityTensor(const MatrixRT& G, const int& phase) const //deprecated
{
return elasticityTensor_(G,phase);
@@ -582,8 +440,6 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
// Getter
///////////////////////////////
Dune::ParameterTree getParameterSet() const {return parameterSet_;}
- // Func2Tensor getElasticityTensor() const {return elasticityTensor_;}
- // Func2TensorParam getElasticityTensor() const {return elasticityTensor_;}
Func2TensorPhase getElasticityTensor() const {return elasticityTensor_;}
MatrixPhase getPrestrainFunction() const {return prestrain_;}
@@ -591,14 +447,16 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
auto getIndicatorFunctionGVF() const {return indicatorFunctionGVF_;} //get as GridViewFunction
auto getIndicatorFunction() const {return indicatorFunction_;}
- //getLameParameters()? .. Throw Exception if isotropic_ = false!
-
///////////////////////////////
// VTK - Writer
///////////////////////////////
- // -----------------------------------------------------------------
- // --- write material (grid)functions to VTK
- // void writeVTKMaterialFunctions(std::array<int, dim> nElements, const int level)
+
+ /**
+ * @brief Write material indicator function as a grid function to VTK.
+ * The Writer uses piecewise constant (P0) interpolation
+ *
+ * @param level GridLevel
+ */
void writeVTKMaterialFunctions(const int level)
{
std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
@@ -622,124 +480,6 @@ Dune::FieldMatrix<double,6,6> setupPhase(const int phase)
return;
};
- // -----------------------------------------------------------------
- // --- write prestrain (grid)functions to VTK
- // void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
- // {
- // std::string outputPath = parameterSet_.get("outputPath", "../../outputs");
- // using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
- // //VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
- // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
-
- // using GridViewVTK = VTKGridType::LeafGridView;
- // const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
-
- // auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
- // auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
-
- // std::vector<double> prestrainFunction_CoeffP0;
- // Functions::interpolate(scalarP0FeBasis, prestrainFunction_CoeffP0, prestrain_);
- // auto prestrainFunction_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, prestrainFunction_CoeffP0);
-
- // std::vector<double> prestrainFunction_CoeffP1;
- // Functions::interpolate(scalarP1FeBasis, prestrainFunction_CoeffP1, prestrain_);
- // auto prestrainFunction_P1 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP1FeBasis, prestrainFunction_CoeffP1);
-
- // VTKWriter<GridView> MaterialVtkWriter(gridView_VTK);
-
- // MaterialVtkWriter.addVertexData(
- // prestrainFunction_P0,
- // VTK::FieldInfo("prestrainFunction_P0", VTK::FieldInfo::Type::scalar, 1));
- // MaterialVtkWriter.addCellData(
- // prestrainFunction_P1,
- // VTK::FieldInfo("prestrainFunction_P1", VTK::FieldInfo::Type::scalar, 1));
-
- // MaterialVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
- // std::cout << "wrote data to file:" + outputPath + "/PrestrainFunctions-level" + std::to_string(level) << std::endl;
- // return;
- // };
-
- // void writeVTKPrestainFunctions(std::array<int, dim> nElements, const int level)
- // {
-// using VTKGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
-// // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{80,80,80});
-// // VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},{40,40,40});
-// VTKGridType grid_VTK({-1.0/2.0, -1.0/2.0, -1.0/2.0},{1.0/2.0, 1.0/2.0, 1.0/2.0},nElements);
-// using GridViewVTK = VTKGridType::LeafGridView;
-// const GridViewVTK gridView_VTK = grid_VTK.leafGridView();
-
-// FTKfillerContainer<dim> VTKFiller;
-// VTKFiller.vtkPrestrainNorm(gridView_VTK, B_Term, "PrestrainBNorm");
-
-// // WORKS Too
-// VTKFiller.vtkProblemCell(gridView_VTK, B_Term, muLocal,"VTKProblemCell");;
-
-
-// // TEST
-// auto scalarP0FeBasis = makeBasis(gridView_VTK,lagrange<0>());
-// auto scalarP1FeBasis = makeBasis(gridView_VTK,lagrange<1>());
-
-// std::vector<double> B_CoeffP0;
-// Functions::interpolate(scalarP0FeBasis, B_CoeffP0, B_Term);
-// auto B_DGBF_P0 = Functions::makeDiscreteGlobalBasisFunction<double>(scalarP0FeBasis, B_CoeffP0);
-
+};
-
-// VTKWriter<GridView> PrestrainVtkWriter(gridView_VTK);
-
-// PrestrainVtkWriter.addCellData(
-// B_DGBF_P0,
-// VTK::FieldInfo("B_P0", VTK::FieldInfo::Type::scalar, 1));
-
-// PrestrainVtkWriter.write(outputPath + "/PrestrainFunctions-level"+ std::to_string(level) );
-// std::cout << "wrote data to file:" + outputPath +"/PrestrainFunctions-level" + std::to_string(level) << std::endl;
-
-// // };
-
- // MatrixRT ElasticityTensorGlobal(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
- // {
- // return MAT(G,x);
- // }
-
- // //--- apply elasticityTensor_ to input Matrix G at position x
- // MatrixRT applyElasticityTensor(const MatrixRT& G, const Domain& x) const //global Version (takes global coordinates)
- // {
- // // Python::Module module = Python::import("material");
- // // auto indicatorFunctionTest = Python::make_function<double>(module.get("indicatorFunction"));
- // // return elasticityTensor_(G,indicatorFunctionTest(x));
- // // auto IFunction = localFunction(indicatorFunction_);
- // // auto IFunction = this->getIndicatorFunction();
- // // auto tmp = Dune::Functions::makeGridViewFunction(indicatorFunction_material_1, gridView_);
- // // auto tmpLocal = LocalF
- // // return elasticityTensor_(G,IFunction(x));
- // // return elasticityTensor_(G,localIndicatorFunction_(x));
- // return elasticityTensor_(G,indicatorFunction_(x));
- // // return elasticityTensor_(G,indicatorFunctionGVF_(x));
- // // int phase = indicatorFunction_(x);
- // // auto tmp = elasticityTensor_(G,phase);
- // // auto tmp = material_1(G,indicatorFunction_(x));
- // // printmatrix(std::cout, material_1(G,indicatorFunction_(x)), "material_1(G,indicatorFunction_(x))", "--");
- // // printmatrix(std::cout, elasticityTensor_(G,phase), "elasticityTensor_(G,phase)", "--");
- // }
-
-
-// template<class Domain>
-// static MatrixRT MAT(const MatrixRT& G, const Domain& x) //unfortunately not possible as a GridViewFunction?
-// {
-// const FieldVector<double,3> mu = {80.0, 80.0, 60.0};
-// const FieldVector<double,3> lambda = {80.0, 80.0, 25.0};
-// double theta=0.25;
-// if (x[0] <(-0.5+theta) and x[2]<(-0.5+theta))
-// return 2.0 * mu[0] * sym(G) + lambda[0] * trace(sym(G)) * Id();
-// else if (x[1]<(-0.5+theta) and x[2]>(0.5-theta))
-// return 2.0 * mu[1] * sym(G) + lambda[1] * trace(sym(G)) * Id();
-// else
-// return 2.0 * mu[2] * sym(G) + lambda[2] * trace(sym(G)) * Id();
-// }
-
-
-
-}; // end class
-
-
-#endif
+#endif
\ No newline at end of file
diff --git a/src/Cell-Problem.cc b/src/Cell-Problem.cc
index 170c46a7..224b0e27 100644
--- a/src/Cell-Problem.cc
+++ b/src/Cell-Problem.cc
@@ -2,7 +2,6 @@
#include <array>
#include <vector>
#include <fstream>
-
#include <fenv.h>
#include <iostream>
@@ -12,13 +11,13 @@
#include <dune/common/parametertreeparser.hh>
#include <dune/common/float_cmp.hh>
#include <dune/common/math.hh>
-
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/yaspgrid.hh>
-// #include <dune/grid/utility/structuredgridfactory.hh> //TEST
#include <dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
#include <dune/istl/matrix.hh>
@@ -30,6 +29,7 @@
#include <dune/istl/spqr.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/io.hh>
+#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/backends/istlvectorbackend.hh>
@@ -42,25 +42,15 @@
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/functions/gridfunctions/gridviewfunction.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
+#include <dune/fufem/dunepython.hh>
-// #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>
#include <dune/microstructure/EffectiveQuantitiesComputer.hh>
#include <dune/microstructure/prestrainedMaterial.hh>
#include <dune/solvers/solvers/umfpacksolver.hh> //TEST
-#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
-
-// #include <dune/fufem/discretizationerror.hh>
-#include <dune/fufem/dunepython.hh>
-// #include <dune/fufem/functions/virtualgridfunction.hh> //TEST
-
-// #include <boost/multiprecision/cpp_dec_float.hpp>
#include <any>
#include <variant>
#include <string>
@@ -112,8 +102,6 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
};
-
-
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
int main(int argc, char *argv[])
@@ -123,7 +111,6 @@ int main(int argc, char *argv[])
Dune::Timer globalTimer;
-
if (argc < 3)
DUNE_THROW(Exception, "Usage: ./Cell-Problem <python path> <python module without extension>");
@@ -144,15 +131,6 @@ int main(int argc, char *argv[])
std::cout << "Input parameters:" << std::endl;
parameterSet.report();
-// ParameterTree parameterSet;
-// if (argc < 2)
-// ParameterTreeParser::readINITree("../../inputs/cellsolver.parset", parameterSet);
-// else
-// {
-// ParameterTreeParser::readINITree(argv[1], parameterSet);
-// ParameterTreeParser::readOptions(argc, argv, parameterSet);
-// }
-
//--- Output setter
std::string baseName = parameterSet.get("baseName", "CellProblem-result");
std::string outputPath = parameterSet.get("outputPath", "../../outputs") ;
@@ -160,27 +138,12 @@ int main(int argc, char *argv[])
std::fstream log;
log.open(outputPath + "/" + baseName + "_log.txt" ,std::ios::out);
-
-// parameterSet.report(log); // short Alternativ
-
- //--- Get Path for Material/Geometry functions
- // std::string geometryFunctionPath = parameterSet.get<std::string>("geometryFunctionPath");
- // //--- Start Python interpreter
- // Python::start();
- // Python::Reference main = Python::import("__main__");
- // Python::run("import math");
- // Python::runStream()
- // << std::endl << "import sys"
- // << std::endl << "sys.path.append('" << geometryFunctionPath << "')"
- // << std::endl;
-
-
constexpr int dim = 3;
// Debug/Print Options
bool print_debug = parameterSet.get<bool>("print_debug", false);
// VTK-write options
-// bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false); //Not implemented yet.
+ // bool write_prestrainFunctions = parameterSet.get<bool>("write_prestrainFunctions", false); //Not implemented yet.
///////////////////////////////////
// Generate the grid
@@ -197,12 +160,10 @@ int main(int argc, char *argv[])
// Create Data Storage
///////////////////////////////////
//--- Storage:: #1 level #2 L2SymError #3 L2SymErrorOrder #4 L2Norm(sym) #5 L2Norm(sym-analytic) #6 L2Norm(phi_1)
-// std::vector<std::variant<std::string, size_t , double>> Storage_Error;
-
+ //std::vector<std::variant<std::string, size_t , double>> Storage_Error;
//--- Storage:: | level | q1 | q2 | q3 | q12 | q13 | q23 | b1 | b2 | b3 |
std::vector<std::variant<std::string, size_t , double>> Storage_Quantities;
-
//--- GridLevel-Loop:
for(size_t level = numLevels[0] ; level <= numLevels[1]; level++)
{
@@ -235,7 +196,7 @@ int main(int argc, char *argv[])
periodicIndices.unifyIndexPair({gridView_CE.indexSet().index(v1)}, {gridView_CE.indexSet().index(v2)});
}
- //--- setup first order periodic Lagrange-Basis
+ //--- Setup first order periodic Lagrange-Basis
auto Basis_CE = makeBasis(
gridView_CE,
power<dim>(
@@ -246,17 +207,15 @@ int main(int argc, char *argv[])
if(print_debug)
std::cout << "power<periodic> basis has " << Basis_CE.dimension() << " degrees of freedom" << std::endl;
-
///////////////////////////////////
// Create prestrained material object
///////////////////////////////////
auto material = prestrainedMaterial(gridView_CE,parameterSet,pyModule);
- // --- get scale ratio
+ // --- Get scale ratio
double gamma = parameterSet.get<double>("gamma",1.0);
std::cout << "scale ratio (gamma) set to : " << gamma << std::endl;
- //------------------------------------------------------------------------------------------------
//--- Compute Correctors
auto correctorComputer = CorrectorComputer(Basis_CE, material, gamma, log, parameterSet);
correctorComputer.assemble();
@@ -278,25 +237,18 @@ int main(int argc, char *argv[])
auto effectiveQuantitiesComputer = EffectiveQuantitiesComputer(correctorComputer,material);
effectiveQuantitiesComputer.computeEffectiveQuantities();
- //--- write material indicator function to VTK
+ //--- Write material indicator function to VTK
if (parameterSet.get<bool>("write_materialFunctions", false))
material.writeVTKMaterialFunctions(level);
- //--- TEST:: Compute Qeff without using the orthogonality (75)...
- // only really makes a difference whenever the orthogonality is not satisfied!
- // std::cout << "----------computeFullQ-----------"<< std::endl; //TEST
- // effectiveQuantitiesComputer.computeFullQ();
-
- //--- get effective quantities
+ //--- Get effective quantities
auto Qeff = effectiveQuantitiesComputer.getQeff();
auto Beff = effectiveQuantitiesComputer.getBeff();
printmatrix(std::cout, Qeff, "Matrix Qeff", "--");
printvector(std::cout, Beff, "Beff", "--");
-
-
- //--- write effective quantities to matlab folder (for symbolic minimization)
+ //--- Write effective quantities to matlab folder (for symbolic minimization)
if(parameterSet.get<bool>("write_toMATLAB", false))
effectiveQuantitiesComputer.writeToMatlab(outputPath);
@@ -306,7 +258,6 @@ int main(int argc, char *argv[])
std::cout<< "q3 : " << std::fixed << Qeff[2][2] << std::endl;
std::cout<< std::fixed << std::setprecision(6) << "q_onetwo=" << Qeff[0][1] << std::endl;
// -------------------------------------------
-
// --- Fill output-Table:
Storage_Quantities.push_back(Qeff[0][0] );
Storage_Quantities.push_back(Qeff[1][1] );
--
GitLab