diff --git a/dune/microstructure/CorrectorComputer.hh b/dune/microstructure/CorrectorComputer.hh
index 576d2afa7963054ec90a7909ad951735c3fb4033..1fbb2ede0009cb59b1d6f229c4e5626ff6834da8 100644
--- a/dune/microstructure/CorrectorComputer.hh
+++ b/dune/microstructure/CorrectorComputer.hh
@@ -1,6 +1,8 @@
#ifndef DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
#define DUNE_MICROSTRUCTURE_CORRECTORCOMPUTER_HH
+#include <map>
+
#include <dune/common/parametertree.hh>
// #include <dune/common/float_cmp.hh>
#include <dune/istl/matrixindexset.hh>
@@ -12,6 +14,8 @@
#include <dune/istl/eigenvalue/test/matrixinfo.hh> // TEST: compute condition Number
#include <dune/solvers/solvers/umfpacksolver.hh>
+#include <dune/microstructure/voigthelper.hh>
+
using namespace MatrixOperations;
using std::shared_ptr;
@@ -45,12 +49,13 @@ protected:
const Material& material_;
+ double gamma_;
+
fstream& log_; // Output-log
const Dune::ParameterTree& parameterSet_;
// const FuncScalar mu_;
// const FuncScalar lambda_;
- double gamma_;
MatrixCT stiffnessMatrix_;
VectorCT load_alpha1_,load_alpha2_,load_alpha3_; //right-hand side(load) vectors
@@ -59,15 +64,14 @@ protected:
VectorCT phi_1_, phi_2_, phi_3_; // Corrector phi_i coefficient vectors
Dune::FieldVector<double,3> m_1_, m_2_, m_3_; // Corrector m_i coefficient vectors
- MatrixRT M1_, M2_, M3_; // (assembled) corrector matrices M_i
- const std::array<MatrixRT*, 3 > mContainer = {&M1_ , &M2_, &M3_};
+ // (assembled) corrector matrices M_i
+ std::array<MatrixRT, 3 > mContainer;
const std::array<VectorCT, 3> phiContainer = {phi_1_,phi_2_,phi_3_};
// ---- Basis for R_sym^{2x2}
- 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 > MatrixBasisContainer_ = {G1_, G2_, G3_};
+ // Could really be constexpr static, but then we would need to write the basis here directly in Voigt notation.
+ // However, I suspect that our Voigt representation may still change in the future.
+ const std::array<VoigtVector<double,3>,3 > matrixBasis_;
Func2Tensor x3G_1_ = [] (const Domain& x) {
return MatrixRT{{1.0*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}; //TODO könnte hier sign übergeben?
@@ -100,6 +104,9 @@ public:
gamma_(gamma),
log_(log),
parameterSet_(parameterSet),
+ matrixBasis_(std::array<VoigtVector<double,3>,3>{matrixToVoigt(Dune::FieldMatrix<double,3,3>({{1, 0, 0}, {0, 0, 0}, {0, 0, 0}})),
+ matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 0, 0}, {0, 1, 0}, {0, 0, 0}})),
+ matrixToVoigt(Dune::FieldMatrix<double,3,3>({{0, 1/std::sqrt(2.0), 0}, {1/std::sqrt(2.0), 0, 0}, {0, 0, 0}}))}),
phiOffset_(basis.size())
{
@@ -153,7 +160,7 @@ public:
// shared_ptr<std::array<VectorRT, 3>> getBasiscontainer(){return make_shared<std::array<VectorRT, 3>>(basisContainer_);}
- auto getMatrixBasiscontainer(){return make_shared<std::array<MatrixRT,3 >>(MatrixBasisContainer_);}
+ 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_;}
@@ -218,7 +225,7 @@ public:
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"
+ assert(arbitraryElementNumber < (std::size_t)basis_.gridView().size(0)); // "arbitraryElementNumber larger than total Number of Elements"
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
@@ -236,7 +243,14 @@ public:
// const localFunction2& lambda
// )
+ /** \brief Compute the stiffness matrix for one element
+ *
+ * \param quadRule The quadrature rule to be used
+ * \param phaseAtQuadPoint The material phase (i.e., the type of material) at each quadrature point
+ */
void computeElementStiffnessMatrix(const typename Basis::LocalView& localView,
+ const Dune::QuadratureRule<double,dim>& quadRule,
+ const std::vector<int>& phaseAtQuadPoint,
ElementMatrixCT& elementMatrix
)
{
@@ -252,14 +266,6 @@ public:
// auto elasticityTensor = material_.getElasticityTensor();
- // auto indicatorFunction = material_.getIndicatorFunction();
- auto localIndicatorFunction = material_.getLocalIndicatorFunction();
- localIndicatorFunction.bind(element);
- // auto indicatorFunctionGVF = material_.getIndicatorFunctionGVF();
- // auto indicatorFunction = localFunction(indicatorFunctionGVF);
- // indicatorFunction.bind(element);
- // Func2int indicatorFunction = *material_.getIndicatorFunction();
-
// LocalBasis-Offset
const int localPhiOffset = localView.size();
@@ -269,28 +275,8 @@ public:
// 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 = Dune::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)
+ for (const auto& quadPoint : quadRule)
{
QPcounter++;
const auto& quadPos = quadPoint.position();
@@ -304,7 +290,7 @@ public:
jacobians[i] = jacobians[i] * geometryJacobianInverse;
// Compute the scaled deformation gradient for each shape function
- std::vector<std::array<MatrixRT,dimworld> > deformationGradient(nSf);
+ std::vector<std::array<VoigtVector<double,dim>,dimworld> > deformationGradient(nSf);
for (size_t i=0; i < nSf; i++)
{
@@ -313,12 +299,12 @@ public:
MatrixRT defGradientV(0);
defGradientV[k] = jacobians[i][0];
- deformationGradient[i][k] = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
+ deformationGradient[i][k] = symVoigt(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
}
}
// Compute the material phase that the current quadrature point is in
- const auto phase = localIndicatorFunction(quadPos);
+ const auto phase = phaseAtQuadPoint[QPcounter-1];
for (size_t l=0; l< dimworld; l++)
for (size_t j=0; j < nSf; j++ )
@@ -371,7 +357,7 @@ public:
// printmatrix(std::cout, material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)), "material_.ElasticityTensor(defGradientU,localIndicatorFunction(quadPos)", "--");
- double energyDensity= scalarProduct(material_.applyElasticityTensor(deformationGradient[i][k],phase),deformationGradient[j][l]);
+ 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));
@@ -412,7 +398,7 @@ public:
- double energyDensityGphi = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),deformationGradient[j][l]);
+ 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;
@@ -461,7 +447,7 @@ public:
- double energyDensityGG = scalarProduct(material_.applyElasticityTensor(basisContainer[m],phase),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]));
@@ -535,14 +521,6 @@ public:
// 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;
@@ -586,19 +564,19 @@ public:
for (size_t k=0; k < dimworld; k++)
{
// Deformation gradient of the ansatz basis function
- MatrixRT defGradientV(0);
- defGradientV[k][0] = jacobians[i][0][0]; // Y
- defGradientV[k][1] = jacobians[i][0][1]; // X2
- // defGradientV[k][2] = (1.0/gamma_)*jacobians[i][0][2]; //
- defGradientV[k][2] = jacobians[i][0][2]; // X3
+ 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
- defGradientV = sym(crossSectionDirectionScaling((1.0/gamma_),defGradientV));
+ 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= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(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));
@@ -624,7 +602,7 @@ public:
// double energyDensityfG= scalarProduct(material_.ElasticityTensorGlobal((-1.0)*forceTerm(quadPos),element.geometry().global(quadPos)),sym(basisContainer[m]));
- double energyDensityfG= scalarProduct(material_.applyElasticityTensor((-1.0)*forceTerm(quadPos),localIndicatorFunction(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]));
@@ -657,6 +635,10 @@ public:
auto localView = basis_.localView();
const int phiOffset = basis_.dimension();
+ // A cache for the element stiffness matrices, if desired
+ bool cacheElementMatrices = true;
+ std::map<std::vector<int>, ElementMatrixCT> elementMatrixCache;
+
// auto muGridF = makeGridViewFunction(mu_, basis_.gridView());
// auto muLocal = localFunction(muGridF);
// auto lambdaGridF = makeGridViewFunction(lambda_, basis_.gridView());
@@ -667,12 +649,51 @@ public:
localView.bind(element);
// muLocal.bind(element);
// lambdaLocal.bind(element);
- const int localPhiOffset = localView.size();
+ const auto localPhiOffset = localView.size();
// Dune::Timer Time;
//std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
+
+ // Determine a suitable quadrature rule
+ const auto& localFiniteElement = localView.tree().child(0).finiteElement();
+ int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
+ const auto& quadRule = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
+
+ // Determine the material phase at each quadrature point
+ auto localIndicatorFunction = material_.getLocalIndicatorFunction();
+ localIndicatorFunction.bind(element);
+
+ std::vector<int> phaseAtQuadPoint(quadRule.size());
+
+ for (std::size_t i=0; i<quadRule.size(); ++i)
+ phaseAtQuadPoint[i] = localIndicatorFunction(quadRule[i].position());
+
ElementMatrixCT elementMatrix;
- // computeElementStiffnessMatrix(localView, elementMatrix, muLocal, lambdaLocal);
- computeElementStiffnessMatrix(localView, elementMatrix);
+
+ // We know that all elements of the grid have the same geometry. Therefore, the element stiffness matrices
+ // depend only on the material phases at the quadrature points. In practice, a few particular phase distributions
+ // dominate (identical phase at all quadrature points). This makes caching the element matrices worthwhile.
+ //
+ // In principle the cache implemented here can get arbitrarily large. If that ever becomes a problem
+ // we need to implement a smarter cache.
+ if (cacheElementMatrices)
+ {
+ auto cachedMatrix = elementMatrixCache.find(phaseAtQuadPoint);
+ if (cachedMatrix != elementMatrixCache.end())
+ {
+ // This copies the matrix. A smarter implementation would avoid this.
+ elementMatrix = (*cachedMatrix).second;
+ }
+ else
+ {
+ computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
+ // This copies the matrix again.
+ elementMatrixCache.insert({phaseAtQuadPoint, elementMatrix});
+ }
+ }
+ else // No caching
+ {
+ computeElementStiffnessMatrix(localView, quadRule, phaseAtQuadPoint, elementMatrix);
+ }
// std::cout << "local assembly time:" << Time.elapsed() << std::endl;
@@ -747,7 +768,7 @@ public:
// lambdaLocal.bind(element);
loadFunctional.bind(element);
- const int localPhiOffset = localView.size();
+ const auto localPhiOffset = localView.size();
// std::cout << "localPhiOffset : " << localPhiOffset << std::endl;
VectorCT elementRhs;
@@ -814,7 +835,7 @@ public:
unsigned int arbitraryLeafIndex = parameterSet_.template get<unsigned int>("arbitraryLeafIndex", 0);
unsigned int arbitraryElementNumber = parameterSet_.template get<unsigned int>("arbitraryElementNumber", 0);
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
- Dune::FieldVector<int,3> row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
+ Dune::FieldVector<std::size_t,3> row = arbitraryComponentwiseIndices(arbitraryElementNumber,arbitraryLeafIndex);
for (int k = 0; k<dim; k++)
{
@@ -902,7 +923,7 @@ public:
auto subtractIntegralMean(VectorCT& coeffVector)
{
- // Substract correct Integral mean from each associated component function
+ // Subtract correct integral mean from each associated component function
auto IM = integralMean(coeffVector);
for(size_t k=0; k<dim; k++)
@@ -1100,30 +1121,27 @@ public:
}
// assemble M_alpha's (actually iota(M_alpha) )
- // MatrixRT M_1(0), M_2(0), M_3(0);
-
- M1_ = 0;
- M2_ = 0;
- M3_ = 0;
+ for (auto& matrix : mContainer)
+ matrix = 0;
for(size_t i=0; i<3; i++)
{
- M1_ += m_1_[i]*MatrixBasisContainer_[i];
- M2_ += m_2_[i]*MatrixBasisContainer_[i];
- M3_ += m_3_[i]*MatrixBasisContainer_[i];
+ mContainer[0] += m_1_[i]*voigtToMatrix(matrixBasis_[i]);
+ mContainer[1] += m_2_[i]*voigtToMatrix(matrixBasis_[i]);
+ mContainer[2] += m_3_[i]*voigtToMatrix(matrixBasis_[i]);
}
std::cout << "--- plot corrector-Matrices M_alpha --- " << std::endl;
- printmatrix(std::cout, M1_, "Corrector-Matrix M_1", "--");
- printmatrix(std::cout, M2_, "Corrector-Matrix M_2", "--");
- printmatrix(std::cout, M3_, "Corrector-Matrix M_3", "--");
+ printmatrix(std::cout, mContainer[0], "Corrector-Matrix M_1", "--");
+ printmatrix(std::cout, mContainer[1], "Corrector-Matrix M_2", "--");
+ printmatrix(std::cout, mContainer[2], "Corrector-Matrix M_3", "--");
log_ << "---------- OUTPUT ----------" << std::endl;
log_ << " --------------------" << std::endl;
- log_ << "Corrector-Matrix M_1: \n" << M1_ << std::endl;
+ log_ << "Corrector-Matrix M_1: \n" << mContainer[0] << std::endl;
log_ << " --------------------" << std::endl;
- log_ << "Corrector-Matrix M_2: \n" << M2_ << std::endl;
+ log_ << "Corrector-Matrix M_2: \n" << mContainer[1] << std::endl;
log_ << " --------------------" << std::endl;
- log_ << "Corrector-Matrix M_3: \n" << M3_ << std::endl;
+ log_ << "Corrector-Matrix M_3: \n" << mContainer[2] << std::endl;
log_ << " --------------------" << std::endl;
@@ -1138,7 +1156,7 @@ public:
}
if(substract_integralMean)
{
- std::cout << " --- substracting integralMean --- " << std::endl;
+ std::cout << " --- subtracting integralMean --- " << std::endl;
subtractIntegralMean(phi_1_);
subtractIntegralMean(phi_2_);
subtractIntegralMean(phi_3_);
@@ -1210,9 +1228,9 @@ public:
computeL2Norm();
computeL2SymGrad();
- 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<< "Frobenius-Norm of M1_: " << mContainer[0].frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M2_: " << mContainer[1].frobenius_norm() << std::endl;
+ std::cout<< "Frobenius-Norm of M3_: " << mContainer[2].frobenius_norm() << std::endl;
}
void computeL2Norm()
@@ -1281,7 +1299,6 @@ public:
auto geometry = element.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- const auto nSf = localFiniteElement.localBasis().size();
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 );
const auto& quad = Dune::QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -1386,21 +1403,17 @@ public:
const double integrationElement = geometry.integrationElement(quadPos);
- auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + *mContainer[b];
+ auto Chi = sym(crossSectionDirectionScaling(1.0/gamma_, DerPhi2(quadPos))) + mContainer[b];
- auto strain1 = DerPhi1(quadPos);
-
-
- strain1 = crossSectionDirectionScaling(1.0/gamma_, strain1);
- strain1 = sym(strain1);
+ const auto strain1 = symVoigt(crossSectionDirectionScaling(1.0/gamma_, DerPhi1(quadPos)));
auto G = matrixFieldG(quadPos);
// auto G = matrixFieldG(e.geometry().global(quadPos)); //TEST
- auto tmp = G + *mContainer[a] + strain1;
+ auto tmp = matrixToVoigt(G + mContainer[a]) + strain1;
- double energyDensity= scalarProduct(material_.applyElasticityTensor(tmp,localIndicatorFunction(quadPos)),sym(Chi));
+ 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);
@@ -1423,7 +1436,7 @@ public:
- // --- Check wheter stiffness matrix is symmetric
+ // --- Check whether stiffness matrix is symmetric
void checkSymmetry()
{
std::cout << " Check Symmetry of stiffness matrix..." << std::endl;
@@ -1433,7 +1446,7 @@ public:
for (const auto& element : elements(basis_.gridView()))
{
localView.bind(element);
- const int localPhiOffset = localView.size();
+ const auto localPhiOffset = localView.size();
for(size_t i=0; i<localPhiOffset; i++)
for(size_t j=0; j<localPhiOffset; j++ )
diff --git a/dune/microstructure/EffectiveQuantitiesComputer.hh b/dune/microstructure/EffectiveQuantitiesComputer.hh
index e7deec3720486f76f10fd09032ef1dbb5c538656..c6e358d4c796396fba5d84d5371d5a975eda3457 100644
--- a/dune/microstructure/EffectiveQuantitiesComputer.hh
+++ b/dune/microstructure/EffectiveQuantitiesComputer.hh
@@ -162,7 +162,7 @@ public:
const auto& quadPos = quadPoint.position();
const double integrationElement = geometry.integrationElement(quadPos);
- auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + *MContainer[a];
+ auto Chi1 = sym(crossSectionDirectionScaling(1.0/gamma, localfun_a(quadPos))) + MContainer[a];
auto G1 = matrixFieldG1(quadPos);
@@ -170,10 +170,10 @@ public:
// auto G1 = matrixFieldG1(e.geometry().global(quadPos)); //TEST
// auto G2 = matrixFieldG2(e.geometry().global(quadPos)); //TEST
- auto X1 = G1 + Chi1;
+ auto X1 = matrixToVoigt(G1 + Chi1);
// auto X2 = G2 + Chi2;
- double energyDensity= scalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(G2));
+ 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() ???
@@ -188,7 +188,7 @@ public:
// printmatrix(std::cout, prestrainPhaseValue_old, "prestrainPhaseValue_old", "--");
// printmatrix(std::cout, prestrainPhaseValue, "prestrainPhaseValue", "--");
// }
- auto value = scalarProduct(material_.applyElasticityTensor(X1,localIndicatorFunction(quadPos)),sym(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;
diff --git a/dune/microstructure/matrix_operations.hh b/dune/microstructure/matrix_operations.hh
index 3eca389fe18d60c9f8d64df98f7aa18b7fb7aa9a..3cd76beddfeda08b1749c784afbe77ee6a9f40bb 100644
--- a/dune/microstructure/matrix_operations.hh
+++ b/dune/microstructure/matrix_operations.hh
@@ -9,26 +9,6 @@ namespace MatrixOperations {
using std::sin;
using std::cos;
- static MatrixRT Id (){
- MatrixRT Id(0);
- for (int i=0;i<3;i++)
- Id[i][i]=1.0;
- return Id;
- }
-
-
- static double norm(VectorRT v){
- return sqrt(pow(v[0],2) + pow(v[1],2) + pow(v[2],2));
- }
-
- static double norm(MatrixRT M){
- return sqrt(
- pow(M[0][0],2) + pow(M[0][1],2) + pow(M[0][2],2)
- + pow(M[1][0],2) + pow(M[1][1],2) + pow(M[1][2],2)
- + pow(M[2][0],2) + pow(M[2][1],2) + pow(M[2][2],2));
- }
-
-
// static MatrixRT sym (MatrixRT M) { // 1/2 (M^T + M)
// MatrixRT ret(0);
// for (int i = 0; i< 3; i++)
@@ -159,10 +139,26 @@ namespace MatrixOperations {
{ R[2][0]*R[2][0], R[2][1]*R[2][1], R[2][2]*R[2][2], R[2][1]*R[2][2], R[2][0]*R[2][2], R[2][0]*R[2][1]},
{2.0*R[1][0]*R[2][0], 2.0*R[1][1]*R[2][1], 2.0*R[1][2]*R[2][2], R[1][1]*R[2][2]+R[1][2]*R[2][1], R[1][0]*R[2][2]+R[1][2]*R[2][0], R[1][0]*R[2][1]+R[1][1]*R[2][0]},
{2.0*R[0][0]*R[2][0], 2.0*R[0][1]*R[2][1], 2.0*R[0][2]*R[2][2], R[0][1]*R[2][2]+R[0][2]*R[2][1], R[0][0]*R[2][2]+R[0][2]*R[2][0], R[0][0]*R[2][1]+R[0][1]*R[2][0]},
- {2.0*R[0][0]*R[0][0], 2.0*R[0][1]*R[1][1], 2.0*R[0][2]*R[1][2], R[0][1]*R[1][2]+R[0][2]*R[1][1], R[0][0]*R[1][2]+R[0][2]*R[1][0], R[0][0]*R[1][1]+R[0][1]*R[1][0]}
+ {2.0*R[0][0]*R[1][0], 2.0*R[0][1]*R[1][1], 2.0*R[0][2]*R[1][2], R[0][1]*R[1][2]+R[0][2]*R[1][1], R[0][0]*R[1][2]+R[0][2]*R[1][0], R[0][0]*R[1][1]+R[0][1]*R[1][0]}
};
}
+
+ /**
+ * @brief Scale last three columns of stiffness matrix by a factor of 2.
+ * Inserting this facotr in the stiffness matrix allows to use the
+ * same Matrix-to-Vector mapping for both the stress and the strain.
+ *
+ * @param S Stiffness matrix
+ */
+ static void scaleStiffnessMatrix(Dune::FieldMatrix<double,6,6>& S){
+ for(size_t i = 0; i<6; i++)
+ for(size_t j = 3; j<6; j++)
+ {
+ S[i][j] = 2.0*S[i][j];
+ }
+ }
+
// static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2){ // ?? Check with Robert
// E1= sym(E1);
// E2 = sym(E2);
@@ -181,7 +177,7 @@ namespace MatrixOperations {
// //
static double linearizedStVenantKirchhoffDensity(double mu, double lambda, MatrixRT E1, MatrixRT E2) // CHANGED
{
- auto t1 = 2.0 * mu * sym(E1) + lambda * trace(sym(E1)) * Id();
+ 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();
@@ -210,7 +206,7 @@ namespace MatrixOperations {
static double QuadraticForm(const double mu, const double lambda, const MatrixRT M){
auto tmp1 = sym(M);
- double tmp2 = norm(tmp1);
+ double tmp2 = tmp1.frobenius_norm();
// double tmp2 = norm(M); //TEST
return lambda*std::pow(trace(M),2) + 2.0*mu*pow( tmp2 ,2);
// return lambda*std::pow(trace(M),2) + 2*mu*pow( norm( sym(M) ),2);
diff --git a/dune/microstructure/prestrainedMaterial.hh b/dune/microstructure/prestrainedMaterial.hh
index 8651815aae24e8822f7cb3e4379cbdfae5c58fe9..deae82234e1553f7c3c0646f9e11b854e3b4fe01 100644
--- a/dune/microstructure/prestrainedMaterial.hh
+++ b/dune/microstructure/prestrainedMaterial.hh
@@ -98,28 +98,16 @@ public:
using MatrixFunc = std::function< MatrixRT(const MatrixRT&) >;
using MatrixPhase = std::function< MatrixRT(const int&) >;
- // using VoigtMatrix = FieldMatrix< double, 6, 6>;
-
protected:
const GridView& gridView_;
const Dune::ParameterTree& parameterSet_;
std::string materialFunctionName_;
- // MatrixFunc L1_;
- // MatrixFunc L2_;
- // MatrixFunc L3_;
-
// Elasticity tensors (in Voigt notation)
- Dune::FieldMatrix<double,6,6> L1_;
- Dune::FieldMatrix<double,6,6> L2_;
- Dune::FieldMatrix<double,6,6> L3_;
- Dune::FieldMatrix<double,6,6> L4_;
+ std::vector<VoigtMatrix<double,3> > L_;
- Func2Tensor prestrain1_;
- Func2Tensor prestrain2_;
- Func2Tensor prestrain3_;
- Func2Tensor prestrain4_;
+ std::vector<Func2Tensor> prestrain_;
Python::Module module_;
@@ -129,7 +117,6 @@ protected:
// Func2Tensor elasticityTensor_;
// Func2TensorParam elasticityTensor_;
Func2TensorPhase elasticityTensor_;
- MatrixPhase prestrain_;
Func2int indicatorFunction_;
@@ -188,9 +175,9 @@ public:
- // L1_ = orthotropic(e1,e2,e3,walnut_parameters);
- // L2_ = orthotropic(e1,e2,e3,Nspruce_parameters);
- // L3_ = orthotropic(e1,e2,e3,Nspruce_parameters);
+ // L_[0] = orthotropic(e1,e2,e3,walnut_parameters);
+ // L_[1] = orthotropic(e1,e2,e3,Nspruce_parameters);
+ // L_[2] = orthotropic(e1,e2,e3,Nspruce_parameters);
setup(materialFunctionName_);
// setup("material");
@@ -200,18 +187,18 @@ public:
localIndicatorFunction_ = localFunction(indicatorFunctionGVF_);
- // L1_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
- // L2_ = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63, 0.37});
- // L3_ = transversely_isotropic(e1,e2,e3, {10.7e3,710 ,500,0.51 ,0.31});
+ // L_[0] = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63 ,0.37});
+ // L_[1] = transversely_isotropic(e1,e2,e3, {11.2e3,1190,960,0.63, 0.37});
+ // L_[2] = transversely_isotropic(e1,e2,e3, {10.7e3,710 ,500,0.51 ,0.31});
- // L1_ = isotropic(mu[0],lambda[0]);
- // L2_ = isotropic(mu[1],lambda[1]);
- // L3_ = isotropic(mu[2],lambda[2]);
+ // L_[0] = isotropic(mu[0],lambda[0]);
+ // L_[1] = isotropic(mu[1],lambda[1]);
+ // L_[2] = isotropic(mu[2],lambda[2]);
- // L1_ = {{lambda[0]+2.0*mu[0], lambda[0], lambda[0], 0.0 , 0.0, 0.0},
+ // L_[0] = {{lambda[0]+2.0*mu[0], lambda[0], lambda[0], 0.0 , 0.0, 0.0},
// {lambda[0], lambda[0]+2*mu[0], lambda[0], 0.0 ,0.0 ,0.0},
// {lambda[0], lambda[0], lambda[0]+2.0*mu[0], 0.0, 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 2.0*mu[0], 0.0, 0.0},
@@ -219,7 +206,7 @@ public:
// { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[0]}
// };
- // L2_ = {{lambda[1]+2.0*mu[1], lambda[1], lambda[1], 0.0 , 0.0, 0.0},
+ // L_[1] = {{lambda[1]+2.0*mu[1], lambda[1], lambda[1], 0.0 , 0.0, 0.0},
// {lambda[1], lambda[1]+2*mu[1], lambda[1], 0.0 ,0.0 ,0.0},
// {lambda[1], lambda[1], lambda[1]+2.0*mu[1], 0.0, 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 2.0*mu[1], 0.0, 0.0},
@@ -227,7 +214,7 @@ public:
// { 0.0 ,0.0, 0.0, 0.0, 0.0, 2.0*mu[1]}
// };
- // L3_ = {{lambda[2]+2.0*mu[2], lambda[2], lambda[2], 0.0 , 0.0, 0.0},
+ // L_[2] = {{lambda[2]+2.0*mu[2], lambda[2], lambda[2], 0.0 , 0.0, 0.0},
// {lambda[2], lambda[2]+2.0*mu[2], lambda[2], 0.0 ,0.0 ,0.0},
// {lambda[2], lambda[2], lambda[2]+2.0*mu[2], 0.0, 0.0, 0.0},
// { 0.0 ,0.0, 0.0, 2.0*mu[2], 0.0, 0.0},
@@ -237,11 +224,11 @@ public:
// printmatrix(std::cout, D_, "D_", "--");
// printmatrix(std::cout, D_inv, "D_inv", "--");
- // printmatrix(std::cout, L1_, "L1_", "--");
- // L1_.leftmultiply(D_inv);
- // printmatrix(std::cout, L1_, "L1_.leftmultiply(D_inv)", "--");
- // L1_.rightmultiply(D_);
- // printmatrix(std::cout, L1_, "L1_.rightmultiply(D_)", "--");
+ // printmatrix(std::cout, L_[0], "L_[0]", "--");
+ // L_[0].leftmultiply(D_inv);
+ // printmatrix(std::cout, L_[0], "L_[0].leftmultiply(D_inv)", "--");
+ // L_[0].rightmultiply(D_);
+ // printmatrix(std::cout, L_[0], "L_[0].rightmultiply(D_)", "--");
// Python::Module module = Python::import(materialFunctionName_);
@@ -251,9 +238,9 @@ public:
// indicatorFunction_ = Python::make_function<int>(module.get("indicatorFunction"));
// indicatorFunction_ = Dune::Functions::makeGridViewFunction(indicatorFunction , gridView_);
// module.get("Phases").toC<int>(Phases_);
- // L1_ = Python::make_function<MatrixRT>(module.get("L1"));
- // L2_ = Python::make_function<MatrixRT>(module.get("L2"));
- // L3_ = Python::make_function<MatrixRT>(module.get("L3"));
+ // L_[0] = Python::make_function<MatrixRT>(module.get("L1"));
+ // L_[1] = Python::make_function<MatrixRT>(module.get("L2"));
+ // L_[2] = Python::make_function<MatrixRT>(module.get("L3"));
// bool isotropic_ = true; // read from module File
// auto Test = Dune::Functions::makeGridViewFunction(MAT<Domain> , gridView_); //not working
}
@@ -278,7 +265,7 @@ Func2Tensor setupPrestrainPhase(Python::Module module, int phase){
module.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
module.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
}
- catch(Dune::Exception)
+ catch(Dune::Exception&)
{
//default frame is used.
}
@@ -342,10 +329,10 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
module.get("phase" + std::to_string(phase) + "_axis").toC<int>(axis);
module.get("phase" + std::to_string(phase) + "_angle").toC<double>(angle);
}
- catch(Dune::Exception)
+ 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 rotaton axis & angle 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;
}
@@ -382,7 +369,7 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
*
* @param name name of the material file being used.
*/
- void setup(const std::string name) // cant use materialFunctionName_ here?
+ 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"));
@@ -391,81 +378,15 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
module.get("Phases").toC<int>(phases);
std::cout << "---- Starting material setup... (Number of Phases: "<< phases << ") " << std::endl;
- switch (phases)
- {
- case 1:
- {
- // std::string phase1_type;
- // module.get("phase1_type").toC<std::string>(phase1_type);
+ L_.resize(phases);
+ prestrain_.resize(phases);
- L1_ = setupPhase(module,1);
- // prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
- prestrain1_ = setupPrestrainPhase(module,1);
- break;
- }
- case 2:
- {
- // std::string phase1_type;
- // module.get("phase1_type").toC<std::string>(phase1_type);
- // std::string phase2_type;
- // module.get("phase2_type").toC<std::string>(phase2_type);
-
- L1_ = setupPhase(module,1);
- L2_ = setupPhase(module,2);
- // prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
- // prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
- prestrain1_ = setupPrestrainPhase(module,1);
- prestrain2_ = setupPrestrainPhase(module,2);
- break;
- }
- case 3:
- {
- // std::string phase1_type;
- // module.get("phase1_type").toC<std::string>(phase1_type);
- // std::string phase2_type;
- // module.get("phase2_type").toC<std::string>(phase2_type);
- // std::string phase3_type;
- // module.get("phase3_type").toC<std::string>(phase3_type);
-
- L1_ = setupPhase(module,1);
- L2_ = setupPhase(module,2);
- L3_ = setupPhase(module,3);
- // prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
- // prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
- // prestrain3_ = Python::make_function<MatrixRT>(module.get("prestrain_phase3"));
- prestrain1_ = setupPrestrainPhase(module,1);
- prestrain2_ = setupPrestrainPhase(module,2);
- prestrain3_ = setupPrestrainPhase(module,3);
- break;
- }
- case 4:
- {
- // std::string phase1_type;
- // module.get("phase1_type").toC<std::string>(phase1_type);
- // std::string phase2_type;
- // module.get("phase2_type").toC<std::string>(phase2_type);
- // std::string phase3_type;
- // module.get("phase3_type").toC<std::string>(phase3_type);
- // std::string phase4_type;
- // module.get("phase4_type").toC<std::string>(phase4_type);
-
- L1_ = setupPhase(module,1);
- L2_ = setupPhase(module,2);
- L3_ = setupPhase(module,3);
- L4_ = setupPhase(module,4);
- // prestrain1_ = Python::make_function<MatrixRT>(module.get("prestrain_phase1"));
- // prestrain2_ = Python::make_function<MatrixRT>(module.get("prestrain_phase2"));
- // prestrain3_ = Python::make_function<MatrixRT>(module.get("prestrain_phase3"));
- // prestrain4_ = Python::make_function<MatrixRT>(module.get("prestrain_phase4"));
- prestrain1_ = setupPrestrainPhase(module,1);
- prestrain2_ = setupPrestrainPhase(module,2);
- prestrain3_ = setupPrestrainPhase(module,3);
- prestrain4_ = setupPrestrainPhase(module,4);
- break;
- }
- default:
- DUNE_THROW(Dune::Exception, " No more than 4 material phases are implemented yet! ");
+ 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;
}
@@ -474,15 +395,15 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
// MATERIAL CLASSES DEFINITIONS:
////////////////////////////////////////////////////////
- //--- Definition of isotropic elasticity tensor (stiffness matrix in voigt notation)
+ /** \brief Isotropic stiffness matrix in Voigt notation but scaled by a factor of 2 in the last three columns. */
Dune::FieldMatrix<double,6,6> isotropic(const double& mu, const double& lambda)
{
- return {{lambda+2.0*mu, lambda , lambda , 0.0 , 0.0 , 0.0 },
- {lambda , lambda+2.0*mu, lambda , 0.0 , 0.0 , 0.0 },
- {lambda , lambda , lambda+2.0*mu, 0.0 , 0.0 , 0.0 },
- { 0.0 , 0.0 , 0.0 , mu , 0.0 , 0.0 },
- { 0.0 , 0.0 , 0.0 , 0.0 , mu , 0.0 },
- { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , mu }
+ return {{lambda+2.0*mu, lambda , lambda , 0.0 , 0.0 , 0.0 },
+ {lambda , lambda+2.0*mu, lambda , 0.0 , 0.0 , 0.0 },
+ {lambda , lambda , lambda+2.0*mu, 0.0 , 0.0 , 0.0 },
+ { 0.0 , 0.0 , 0.0 , 2.0*mu, 0.0 , 0.0 },
+ { 0.0 , 0.0 , 0.0 , 0.0 , 2.0*mu, 0.0 },
+ { 0.0 , 0.0 , 0.0 , 0.0 , 0.0 , 2.0*mu }
};
}
//-------------------------------------------------------------------------------
@@ -555,12 +476,14 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
(Compliance.leftmultiply(C_rot)).rightmultiply(C_rot.transposed());
}
-
printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
- // invert to obtain elasticity/stiffness-Tensor
+ // invert Compliance matrix to obtain elasticity/stiffness matrix
Compliance.invert();
// printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
+ scaleStiffnessMatrix(Compliance);
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor scaled", "--");
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
return Compliance;
}
//-------------------------------------------------------------------------------
@@ -627,7 +550,10 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
//--- return elasticity tensor
Compliance.invert();
- // printmatrix(std::cout, C, "C after .invert()", "--");
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
+
+ scaleStiffnessMatrix(Compliance);
+
return Compliance;
}
//-------------------------------------------------------------------------------
@@ -667,85 +593,35 @@ Dune::FieldMatrix<double,6,6> setupPhase(Python::Module module, int phase)
}
printmatrix(std::cout, Compliance, "Compliance matrix used:", "--");
- //--- return elasticity tensor (in Voigt notation)
+
+
+
+ //--- return stiffness/elasticity tensor (in Voigt notation)
Compliance.invert();
- // printmatrix(std::cout, C, "C after .invert()", "--");
+ // printmatrix(std::cout, Compliance, "Stiffness Tensor", "--");
+
+ scaleStiffnessMatrix(Compliance);
+
return Compliance;
}
//-------------------------------------------------------------------------------
- MatrixRT applyElasticityTensor(const MatrixRT& G, const int& phase) const
+ VoigtVector<double,3> applyElasticityTensor(const VoigtVector<double,3>& G, const int& phase) const
{
-
- VoigtVector<double,3> G_tmp = strainToVoigt(G);
VoigtVector<double,3> out(0);
- switch (phase)
- {
- case 1:
- {
- // printmatrix(std::cout, L1_, "L1_", "--");
- L1_.mv(G_tmp,out);
- break;
- }
- case 2:
- {
- // printmatrix(std::cout, L2_, "L2_", "--");
- L2_.mv(G_tmp,out);
- break;
- }
- case 3:
- {
- // printmatrix(std::cout, L3_, "L3_", "--");
- L3_.mv(G_tmp,out);
- break;
- }
- case 4:
- {
- // printmatrix(std::cout, L4_, "L4_", "--");
- L4_.mv(G_tmp,out);
- break;
- }
- default:
- DUNE_THROW(Dune::Exception, "Phase number not implemented.");
- break;
- }
-
- return voigtToStress(out);
+ // printmatrix(std::cout, L_[0], "L_[0]", "--");
+ L_[phase-1].mv(G,out);
+ return out;
}
////////////////////////////////////////////////////////////////////////////////////////
MatrixRT prestrain(const int& phase,const Domain& x) const
{
- switch (phase)
- {
- case 1:
- {
- return prestrain1_(x);
- break;
- }
- case 2:
- {
- return prestrain2_(x);
- break;
- }
- case 3:
- {
- return prestrain3_(x);
- break;
- }
- case 4:
- {
- return prestrain3_(x);
- break;
- }
- default:
- DUNE_THROW(Dune::Exception, "Phase number not implemented.");
- break;
- }
+ return prestrain_[phase-1](x);
}
diff --git a/dune/microstructure/voigthelper.hh b/dune/microstructure/voigthelper.hh
index db0bf17c49e3049950de77900f77c5d16ee826d0..5058f710ec2b21a535b8da7f51aec7082d498646 100644
--- a/dune/microstructure/voigthelper.hh
+++ b/dune/microstructure/voigthelper.hh
@@ -16,28 +16,49 @@ template<typename T, int dim>
using VoigtMatrix = Dune::FieldMatrix<T,dim*(dim+1)/2,dim*(dim+1)/2>;
+/**
+ * @brief Voigt-notation(https://de.wikipedia.org/wiki/Voigtsche_Notation) distinguishes
+ * in the transformation from Matrix to Vector between stresses and strains. The transformation
+ * for strains features an additional factor 2 for the non-diagonal entries. In order to avoid
+ * the use of different data structures for both stresses & strains we use the same Matrix-to-Vector
+ * mapping ('matrixToVoigt') and incorporate the factors in suitable places. namely:
+ * - The Stiffness matrix of the constitutive relation get scaled by a factor of 2 in the last three columns
+ * - The 'voigtScalarProduct' scales the lase three products by a factor of 2
+ */
+
template<typename T>
-static VoigtVector<T,3> strainToVoigt(const Dune::FieldMatrix<T,3,3>& G)
+VoigtVector<T,3> matrixToVoigt(const Dune::FieldMatrix<T,3,3>& matrix)
{
- return {G[0][0], G[1][1], G[2][2], 2.0*G[1][2], 2.0*G[0][2], 2.0*G[0][1]};
+ return {matrix[0][0], matrix[1][1], matrix[2][2], matrix[1][2], matrix[0][2], matrix[0][1]};
}
template<typename T>
-static Dune::FieldMatrix<T,3,3> voigtToStrain(const VoigtVector<T,3>& v)
+Dune::FieldMatrix<T,3,3> voigtToMatrix(const VoigtVector<T,3>& v)
{
- return {{v[0], v[5]/2.0, v[4]/2.0}, {v[5]/2.0, v[1], v[3]/2.0}, {v[4]/2.0, v[3]/2.0, v[2]}};
+ return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
}
+/** \brief Interpret Voigt vectors as symmetric matrices and compute their Frobenius scalar voigtScalarProduct
+ */
template<typename T>
-static VoigtVector<T,3> stressToVoigt(const Dune::FieldMatrix<T,3,3>& G)
+T voigtScalarProduct(const VoigtVector<T,3>& a, const VoigtVector<T,3>& b)
{
- return {G[0][0], G[1][1], G[2][2], G[1][2], G[0][2], G[0][1]};
+ return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] + 2.0*a[3]*b[3] + 2.0*a[4]*b[4] + 2.0*a[5]*b[5];
}
+/** \brief Return the symmetric part of a matrix, in Voigt notation */
template<typename T>
-static Dune::FieldMatrix<T,3,3> voigtToStress(const VoigtVector<T,3>& v)
+VoigtVector<T,3> symVoigt(const Dune::FieldMatrix<T,3,3>& matrix)
{
- return {{v[0], v[5], v[4]}, {v[5], v[1], v[3]}, {v[4], v[3], v[2]}};
+ VoigtVector<T,3> result = {matrix[0][0],
+ matrix[1][1],
+ matrix[2][2],
+ // The 0.5 is part of the sym operator, the 2.0 is from the conversion to Voigt notation
+ 0.5 * (matrix[1][2] + matrix[2][1]),
+ 0.5 * (matrix[0][2] + matrix[2][0]),
+ 0.5 * (matrix[0][1] + matrix[1][0])};
+
+ return result;
}
#endif // DUNE_MICROSTRUCTURE_VOIGTHELPER_HH