diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index d9b1d490bb1173b15f5b9338caefe833e7cb05af..0fb19495532575032af50a86c03fa3881117d8a8 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -16,7 +16,8 @@ cellDomain = 1
#nElements_Cell = 20 20 20 # number elements in each direction (y1 y2 x3)
#nElements_Cell = 100 100 2
-nElements_Cell = 4 4 4
+nElements_Cell = 10 10 10
+#nElements_Cell = 4 4 4
width = 1.0
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index d188d12287208c1c3166454dba5109697b5d243d..9fee3c166847f182f2e0b991299693e27f6e298b 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -61,7 +61,7 @@ auto arbitraryComponentwiseIndices(const Basis& basis,
const int leafIdx
)
{
- // (Local Approach -- works for non Lagrangian-Basis too)
+ // (Local Approach -- works for non Lagrangian-Basis too)
// Input : elementNumber & localIdx
// Output : determine all Component-Indices that correspond to that basis-function
auto localView = basis.localView();
@@ -74,7 +74,7 @@ auto arbitraryComponentwiseIndices(const Basis& basis,
if(elementCounter == elementNumber) // get arbitraryIndex(global) for each Component ..alternativ:gridView.indexSet
{
localView.bind(element);
-
+
for (int k = 0; k < 3; k++)
{
auto rowLocal = localView.tree().child(k).localIndex(leafIdx); //Input: LeafIndex! TODO braucht hier (Inverse ) Local-to-Leaf Idx Map
@@ -127,8 +127,8 @@ void computeElementStiffnessMatrix(const LocalView& localView,
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
-// std::cout << "localView.size(): " << localView.size() << std::endl;
-// std::cout << "nSf : " << nSf << std::endl;
+// 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...
@@ -143,7 +143,7 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// printmatrix(std::cout, basisContainer[2] , "G_3", "--");
////////////////////////////////////////////////////
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -295,14 +295,14 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// for (size_t k=0; k < nCompo; k++)
// for (size_t i=0; i < nSf; i++)
// {
-//
+//
// // (scaled) Deformation gradient of the ansatz basis function
// MatrixRT defGradientU(0);
// defGradientU[k][0] = gradients[i][0]; // Y
// defGradientU[k][1] = gradients[i][1]; //X2
// defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
-//
-//
+//
+//
// // symmetric Gradient (Elastic Strains)
// MatrixRT strainU;
// for (int ii=0; ii<nCompo; ii++)
@@ -310,34 +310,34 @@ void computeElementStiffnessMatrix(const LocalView& localView,
// {
// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
// }
-//
+//
// // St.Venant-Kirchhoff stress
// MatrixRT stressU(0);
// stressU.axpy(2*mu(quadPos), strainU); //this += 2mu *strainU
-//
+//
// double trace = 0;
// for (int ii=0; ii<nCompo; ii++)
// trace += strainU[ii][ii];
-//
+//
// for (int ii=0; ii<nCompo; ii++)
// stressU[ii][ii] += lambda(quadPos) * trace;
-//
+//
// for( size_t n=0; n<3; n++)
// {
-//
+//
// // < L sym[D_gamma*nabla phi_i], G_j >
// double energyDensityPhiG = 0;
// for (int ii=0; ii<nCompo; ii++)
// for (int jj=0; jj<nCompo; jj++)
// energyDensityPhiG += stressU[ii][jj] * basisContainer[n][ii][jj]; // "phi*m"-part
-//
+//
// size_t col = localView.tree().child(k).localIndex(i);
-//
+//
// elementMatrix[localPhiOffset+n][col] += energyDensityPhiG * quadPoint.weight() * integrationElement;
-//
+//
// }
// }
-//
+//
@@ -382,12 +382,12 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
// OccupationPattern:
// | phi*phi m*phi |
// | phi *m m*m |
-
+
auto localView = basis.localView();
const int phiOffset = basis.dimension();
- nb.resize(basis.size()+3, basis.size()+3);
+ nb.resize(basis.size()+3, basis.size()+3);
for (const auto& element : elements(basis.gridView()))
{
@@ -421,31 +421,31 @@ void getOccupationPattern(const Basis& basis, MatrixIndexSet& nb, ParameterSet&
}
}
//////////////////////////////////////////////////////////////////
- // setOneBaseFunctionToZero
+ // setOneBaseFunctionToZero
//////////////////////////////////////////////////////////////////
- if(parameterSet.template get<bool>("set_IntegralZero", true)){
+ if(parameterSet.template get<bool>("set_IntegralZero", true)){
FieldVector<int,3> row;
- unsigned int arbitraryLeafIndex = parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
+ unsigned int arbitraryLeafIndex = parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
unsigned int arbitraryElementNumber = parameterSet.template get<unsigned int>("arbitraryElementNumber", 0);
-
-
-
-// assert(arbitraryLeafIndex < localView.size() ); // "localIndex is larger than #shapeFunctions" TODO ERROR
-
-
+
+
+
+// assert(arbitraryLeafIndex < localView.size() ); // "localIndex is larger than #shapeFunctions" TODO ERROR
+
+
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 );
-
+ const auto nSf = localFiniteElement.localBasis().size();
+ assert(arbitraryLeafIndex < nSf );
+
std::cout << "arbitraryElementNumber:" << arbitraryElementNumber << std::endl;
std::cout << "basis.gridView().size(0:" << basis.gridView().size(0) << std::endl;
assert(arbitraryElementNumber < basis.gridView().size(0)); // "arbitraryElementNumber larger than total Number of Elements"
-
+
//Determine 3 global indices (one for each component of an arbitrary local FE-function)
row = arbitraryComponentwiseIndices(basis,arbitraryElementNumber,arbitraryLeafIndex);
-
+
printvector(std::cout, row, "row" , "--");
-
+
for (int k = 0; k<3; k++)
nb.add(row[k],row[k]);
}
@@ -499,12 +499,12 @@ void computeElementLoadVector( const LocalView& localView,
MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 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)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1); // für simplex
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
- for (const auto& quadPoint : quad)
+ for (const auto& quadPoint : quad)
{
@@ -553,7 +553,7 @@ void computeElementLoadVector( const LocalView& localView,
double energyDensity = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- energyDensity += stressV[ii][jj] *forceTerm(quadPos)[ii][jj];
+ energyDensity += stressV[ii][jj] *forceTerm(quadPos)[ii][jj];
size_t row = localView.tree().child(k).localIndex(i);
@@ -579,7 +579,7 @@ void computeElementLoadVector( const LocalView& localView,
double energyDensityfG = 0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
+ energyDensityfG += stressG[ii][jj] * forceTerm(quadPos)[ii][jj];
elementRhs[localPhiOffset+m] += energyDensityfG * quadPoint.weight() * integrationElement;
}
@@ -719,8 +719,8 @@ void assembleCellLoad(const Basis& basis,
{
b[phiOffset+m] += elementRhs[localPhiOffset+m];
}
-
-
+
+
// printvector(std::cout, b, "b", "--");
}
@@ -753,10 +753,10 @@ auto energy(const Basis& basis,
using GridView = typename Basis::GridView;
using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
-
-
-
- // TEST
+
+
+
+ // TEST
// FieldVector<double,3> testvector = {1.0 , 1.0 , 1.0};
// printmatrix(std::cout, matrixFieldFuncB(testvector) , "matrixFieldB(testvector) ", "--");
@@ -776,14 +776,14 @@ auto energy(const Basis& basis,
auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+3; // TEST
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+3; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(e.type(), orderQR);
for (const auto& quadPoint : quad) {
// for (size_t pt=0; pt < quad.size(); pt++) {
-
-
+
+
const FieldVector<double,dim>& quadPos = quadPoint.position();
// const Domain& quadPos = quad[pt].position();
const double integrationElement = geometry.integrationElement(quadPos);
@@ -836,7 +836,7 @@ void setOneBaseFunctionToZero(const Basis& basis,
constexpr int dim = 3;
FieldVector<int,3> row;
- unsigned int arbitraryLeafIndex = parameterSet.template get<unsigned int>("arbitraryLeafIndex", 0);
+ 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)
@@ -844,7 +844,7 @@ void setOneBaseFunctionToZero(const Basis& basis,
for (int k = 0; k < dim; k++)
{
- load_alpha1[row[k]] = 0.0;
+ load_alpha1[row[k]] = 0.0;
load_alpha2[row[k]] = 0.0;
load_alpha3[row[k]] = 0.0;
auto cIt = stiffnessMatrix[row[k]].begin();
@@ -865,8 +865,8 @@ auto childToIndexMap(const Basis& basis,
const int k
)
{
- // Input : child/component
- // Output : determine all Indices that correspond to that component
+ // Input : child/component
+ // Output : determine all Indices that correspond to that component
auto localView = basis.localView();
@@ -880,11 +880,11 @@ auto childToIndexMap(const Basis& basis,
{
localView.bind(element);
const auto& localFiniteElement = localView.tree().child(k).finiteElement(); // macht keinen Unterschied ob hier k oder 0..
- const auto nSf = localFiniteElement.localBasis().size(); //
+ const auto nSf = localFiniteElement.localBasis().size(); //
for(size_t j=0; j<nSf; j++)
{
- auto Localidx = localView.tree().child(k).localIndex(j); // local indices
+ auto Localidx = localView.tree().child(k).localIndex(j); // local indices
r.push_back(localView.index(Localidx)); // global indices
}
}
@@ -912,8 +912,8 @@ auto integralMean(const Basis& basis,
)
{
// computes integral mean of given LocalFunction
- constexpr int dim = 3;
-
+ constexpr int dim = 3;
+
auto GVFunction = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
auto localfun = localFunction(GVFunction);
@@ -928,7 +928,7 @@ auto integralMean(const Basis& basis,
localView.bind(element);
localfun.bind(element);
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
-
+
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const QuadratureRule<double, dim>& quad = QuadratureRules<double, dim>::rule(element.type(), orderQR);
@@ -936,7 +936,7 @@ auto integralMean(const Basis& basis,
{
const double integrationElement = element.geometry().integrationElement(quadPoint.position());
area += quadPoint.weight() * integrationElement;
-
+
const auto& quadPos = quadPoint.position();
r += localfun(quadPos) * quadPoint.weight() * integrationElement;
}
@@ -995,12 +995,12 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
-
-
-
-
-
- /*
+
+
+
+
+
+ /*
template<class Basis, class Vector, class MatrixFunction, class Domain> //TODO
double evalSymGrad(const Basis& basis,
@@ -1009,7 +1009,7 @@ double evalSymGrad(const Basis& basis,
Domain& x // evaluation Point in reference coordinates
)
{
-
+
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -1023,21 +1023,21 @@ double evalSymGrad(const Basis& basis,
-
-
+
+
for (const auto& element : elements(basis.gridView()))
{
-
+
localView.bind(element);
auto geometry = element.geometry();
-
+
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
-
+
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(x);
-
+
std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
localFiniteElement.localBasis().evaluateJacobian(x,
@@ -1052,22 +1052,22 @@ double evalSymGrad(const Basis& basis,
MatrixRT defGradientU(0);
-
-
+
+
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++) //Error: these are local fcts!
{
-
- size_t localIdx = localView.tree().child(k).localIndex(i);
+
+ size_t localIdx = localView.tree().child(k).localIndex(i);
size_t globalIdx = localView.index(localIdx);
// (scaled) Deformation gradient of the ansatz basis function
defGradientU[k][0] += coeffVector[globalIdx]* gradients[i][0]; // Y
defGradientU[k][1] += coeffVector[globalIdx]* gradients[i][1]; //X2
defGradientU[k][2] += coeffVector[globalIdx]*(1.0/gamma)*gradients[i][2]; //X3
-
-
-
+
+
+
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
for (int ii=0; ii<nCompo; ii++)
@@ -1075,8 +1075,8 @@ double evalSymGrad(const Basis& basis,
{
strainU[ii][jj] = coeffVector[globalIdx] * 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
}
-
-
+
+
}
}
@@ -1093,7 +1093,7 @@ double computeL2SymError(const Basis& basis,
const MatrixFunction& matrixFieldFunc
)
{
- auto error = 0.0;
+ double error = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -1111,9 +1111,10 @@ double computeL2SymError(const Basis& basis,
localView.bind(element);
matrixField.bind(element);
auto geometry = element.geometry();
-
+
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
+// std::cout << "print nSf:" << nSf << std::endl;
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
@@ -1138,84 +1139,83 @@ double computeL2SymError(const Basis& basis,
MatrixRT defGradientU(0), defGradientV(0);
-
-
+
+
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
-
{
- for (size_t l=0; l< nCompo; l++)
+ for (size_t l=0; l < nCompo; l++)
for (size_t j=0; j < nSf; j++ )
{
-
- size_t localIdx1 = localView.tree().child(k).localIndex(i);
+
+ size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
size_t globalIdx1 = localView.index(localIdx1);
- size_t localIdx2 = localView.tree().child(l).localIndex(j);
+ size_t localIdx2 = localView.tree().child(l).localIndex(j);
size_t globalIdx2 = localView.index(localIdx2);
// (scaled) Deformation gradient of the ansatz basis function
- defGradientU[k][0] = gradients[i][0]; // Y
+ defGradientU[k][0] = gradients[i][0]; // Y //hier i:leaf oder localIdx?
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
-
+
defGradientV[l][0] = gradients[j][0]; // Y
defGradientV[l][1] = gradients[j][1]; //X2
defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
-
+
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0), strainV(0);
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
- strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
- strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
+ strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
+ strainV[ii][jj] = 0.5 * (defGradientV[ii][jj] + defGradientV[jj][ii]);
}
// Local energy density: stress times strain
double tmp = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
+ tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
error += tmp * quadPoint.weight() * integrationElement;
}
}
-
+
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
- size_t localIdx1 = localView.tree().child(k).localIndex(i);
+ size_t localIdx1 = localView.tree().child(k).localIndex(i);
size_t globalIdx1 = localView.index(localIdx1);
// (scaled) Deformation gradient of the ansatz basis function
defGradientU[k][0] = gradients[i][0]; // Y
defGradientU[k][1] = gradients[i][1]; //X2
defGradientU[k][2] = (1.0/gamma)*gradients[i][2]; //X3
-
+
// symmetric Gradient (Elastic Strains)
MatrixRT strainU(0);
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
- strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
+ strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
double tmp = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- tmp += (-2.0)*coeffVector[globalIdx1]*strainU[ii][jj] * matrixField(quadPos)[ii][jj];
-
+ tmp += (-2.0)*coeffVector[globalIdx1]*strainU[ii][jj] * matrixField(quadPos)[ii][jj];
+
error += tmp * quadPoint.weight() * integrationElement;
}
-
+
double tmp = 0.0;
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- tmp += matrixField(quadPos)[ii][jj] * matrixField(quadPos)[ii][jj];
-
+ tmp += matrixField(quadPos)[ii][jj] * matrixField(quadPos)[ii][jj];
+
error += tmp * quadPoint.weight() * integrationElement;
-
+
}
}
return sqrt(error);
@@ -1229,7 +1229,7 @@ double computeL2SymError(const Basis& basis,
-/*
+
template<class Basis, class Vector, class MatrixFunction>
double computeL2ErrorOld(const Basis& basis,
Vector& coeffVector,
@@ -1237,7 +1237,7 @@ double computeL2ErrorOld(const Basis& basis,
const MatrixFunction& matrixFieldFunc
)
{
-
+
auto error = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -1254,19 +1254,19 @@ double computeL2ErrorOld(const Basis& basis,
for (const auto& element : elements(basis.gridView()))
{
-
+
localView.bind(element);
matrixField.bind(element);
auto geometry = element.geometry();
-
+
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); // Unterscheidung children notwendig?
const auto nSf = localFiniteElement.localBasis().size();
-
-// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
+
+// int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1)+5; // TEST
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
+
for (const auto& quadPoint : quad)
{
const auto& quadPos = quadPoint.position();
@@ -1286,13 +1286,13 @@ double computeL2ErrorOld(const Basis& basis,
MatrixRT defGradientU(0);
-
-
+
+
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
-
- size_t localIdx = localView.tree().child(k).localIndex(i);
+
+ size_t localIdx = localView.tree().child(k).localIndex(i);
size_t globalIdx = localView.index(localIdx);
// (scaled) Deformation gradient of the ansatz basis function
@@ -1306,7 +1306,7 @@ double computeL2ErrorOld(const Basis& basis,
for (int jj=0; jj<nCompo; jj++)
{
strainU[ii][jj] = coeffVector[globalIdx] * 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
-// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient //TEST
+// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient //TEST
}
// Local energy density: stress times strain
@@ -1315,27 +1315,27 @@ double computeL2ErrorOld(const Basis& basis,
for (int jj=0; jj<nCompo; jj++)
tmp+= std::pow(strainU[ii][jj] - matrixField(quadPos)[ii][jj],2);
// tmp+= std::pow((coeffVector[globalIdx]*strainU[ii][jj]) - matrixField(quadPos)[ii][jj] ,2);
-
+
error += tmp * quadPoint.weight() * integrationElement;
}
}
}
return sqrt(error);
-}*/
+}
/*
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();
@@ -1370,12 +1370,12 @@ for (const auto& e : elements(basis_.gridView()))
MatrixRT defGradient2(0);
defGradient2[l] = gradients[j];
- size_t localIdx1 = localView.tree().child(k).localIndex(i);
+ 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 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++)
@@ -1383,7 +1383,7 @@ for (const auto& e : elements(basis_.gridView()))
elementScp += scp;
}
-
+
ret += elementScp * quadPoint.weight() * integrationElement;
} //qp
} //e
@@ -1414,7 +1414,7 @@ int main(int argc, char *argv[])
ParameterTreeParser::readINITree(argv[1], parameterSet);
ParameterTreeParser::readOptions(argc, argv, parameterSet);
}
-
+
// output setter
std::string outputPath = parameterSet.get("outputPath", "../../outputs/output.txt");
std::fstream log;
@@ -1424,7 +1424,7 @@ int main(int argc, char *argv[])
debugLog.open("../../outputs/debugLog.txt",std::ios::out);
// log << "Writing this to a file. \n";
// log.close();
-
+
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -1443,7 +1443,7 @@ int main(int argc, char *argv[])
///////////////////////////////////
- // Get Prestrain Parameters
+ // Get Prestrain Parameters
///////////////////////////////////
unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", 2);
double rho1 = parameterSet.get<double>("rho1", 1.0);
@@ -1456,13 +1456,13 @@ int main(int argc, char *argv[])
log << "prestrain imp: " << prestraintype << "\nrho1 = " << rho1 << "\nrho2 = " << rho2 << std::endl;
log << "alpha: " << alpha << std::endl;
- log << "gamma: " << gamma << std::endl;
-
+ log << "gamma: " << gamma << std::endl;
+
///////////////////////////////////
// Generate the grid
///////////////////////////////////
-
+
//Corrector Problem Domain:
unsigned int cellDomain = parameterSet.get<unsigned int>("cellDomain", 1);
// (shifted) Domain (-1/2,1/2)^3
@@ -1480,7 +1480,7 @@ int main(int argc, char *argv[])
log << "Number of Elements in each direction: " << nElements << std::endl;
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
- CellGridType grid_CE(lower,upper,nElements);
+ CellGridType grid_CE(lower,upper,nElements);
using GridView = CellGridType::LeafGridView;
const GridView gridView_CE = grid_CE.leafGridView();
@@ -1503,7 +1503,7 @@ int main(int argc, char *argv[])
double mu2 = beta*mu1;
double lambda1 = parameterSet.get<double>("lambda1",0.0);;
double lambda2 = beta*lambda1;
-
+
///////////////////////////////////
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
@@ -1524,22 +1524,22 @@ int main(int argc, char *argv[])
return lambda2;
};
- auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
+ auto muGridF = Dune::Functions::makeGridViewFunction(muTerm, gridView_CE);
auto muLocal = localFunction(muGridF);
auto lambdaGridF = Dune::Functions::makeGridViewFunction(lambdaTerm, gridView_CE);
auto lambdaLocal = localFunction(lambdaGridF);
-
+
log << "beta: " << beta << std::endl;
- log << "material parameters: " << std::endl;
+ log << "material parameters: " << std::endl;
log << "mu1: " << mu1 << "\nmu2: " << mu2 << std::endl;
log << "lambda1: " << lambda1 <<"\nlambda2:" << lambda2 << std::endl;
-
+
///////////////////////////////////
// --- Choose Solver ---
// 1 : CG-Solver
// 2 : GMRES
- // 3 : QR
+ // 3 : QR
///////////////////////////////////
unsigned int Solvertype = parameterSet.get<unsigned int>("Solvertype", 1);
bool write_corrector_phi1 = parameterSet.get<bool>("write_corrector_phi1", false);
@@ -1547,8 +1547,8 @@ int main(int argc, char *argv[])
bool write_corrector_phi3 = parameterSet.get<bool>("write_corrector_phi3", false);
bool write_L2Error = parameterSet.get<bool>("write_L2Error", false);
bool write_IntegralMean = parameterSet.get<bool>("write_IntegralMean", false);
-
-
+
+
/////////////////////////////////////////////////////////
// Choose a finite element space for Cell Problem
/////////////////////////////////////////////////////////
@@ -1556,12 +1556,12 @@ int main(int argc, char *argv[])
Functions::BasisFactory::Experimental::PeriodicIndexSet periodicIndices;
- int equivCounter = 0; // TODO remove
+ int equivCounter = 0; // TODO remove
// Don't do the following in real life: It has quadratic run-time in the number of vertices.
for (const auto& v1 : vertices(gridView_CE))
{
-// std::cout << " ---------------------------------------" << std::endl;
+// std::cout << " ---------------------------------------" << std::endl;
for (const auto& v2 : vertices(gridView_CE))
if (equivalent(v1.geometry().corner(0), v2.geometry().corner(0)))
{
@@ -1572,10 +1572,10 @@ int main(int argc, char *argv[])
// std::cout << "equivCounter:" << equivCounter << std::endl;
// std::cout << "Number of Elements in each direction: " << nE << std::endl;
// std::cout << "Number ofNodes: " << gridView_CE.size(dim) << std::endl;
-
-
- auto perTest = periodicIndices.indexPairSet(); // TODO remove
-
+
+
+ auto perTest = periodicIndices.indexPairSet(); // TODO remove
+
auto Basis_CE = makeBasis(
gridView_CE,
power<dim>( //lagrange<1>(),
@@ -1587,10 +1587,10 @@ int main(int argc, char *argv[])
std::cout << "size feBasis: " << Basis_CE.size() << std::endl;
std::cout << "basis.dimension()" << Basis_CE.dimension() <<std::endl;
log << "size of FiniteElementBasis: " << Basis_CE.size() << std::endl;
-
+
const int phiOffset = Basis_CE.size();
debugLog << "phiOffset: "<< phiOffset << std::endl;
-
+
/////////////////////////////////////////////////////////
// Data structures: Stiffness matrix and right hand side vector
/////////////////////////////////////////////////////////
@@ -1620,7 +1620,7 @@ int main(int argc, char *argv[])
};
Func2Tensor x3G_3 = [] (const Domain& x) {
- return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
Func2Tensor x3G_1neg = [x3G_1] (const Domain& x) {
@@ -1634,29 +1634,29 @@ int main(int argc, char *argv[])
Func2Tensor x3G_3neg = [x3G_3] (const Domain& x) {
return -1.0*x3G_3(x);
};
-
-
+
+
///////////////////////////////////////////////
// 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}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
// printmatrix(std::cout, basisContainer[2] , "´G_3", "--");
// log << basisContainer[0] << std::endl;
-
+
//////////////////////////////////////////////////////////////////////
- // Determine global indices of components for arbitrary (local) index
+ // Determine global indices of components for arbitrary (local) index
//////////////////////////////////////////////////////////////////////
- int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
+ int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
- FieldVector<int,3> row;
+ FieldVector<int,3> row;
row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
printvector(std::cout, row, "row" , "--");
@@ -1667,22 +1667,22 @@ int main(int argc, char *argv[])
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha1 ,x3G_1neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha2 ,x3G_2neg);
assembleCellLoad(Basis_CE, muLocal, lambdaLocal,gamma, load_alpha3 ,x3G_3neg);
-
+
// printmatrix(std::cout, stiffnessMatrix_CE, "StiffnessMatrix", "--");
// printvector(std::cout, load_alpha1, "load_alpha1", "--");
-
-
+
+
//////////////////////////////////////////////////////////////////////
- // Determine global indices of components for arbitrary (local) index
+ // Determine global indices of components for arbitrary (local) index
//////////////////////////////////////////////////////////////////////
-// int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
+// int arbitraryLeafIndex = parameterSet.get<unsigned int>("arbitraryLeafIndex", 0); // localIdx..assert Number < #ShapeFcts .. also input Element?
// int arbitraryElementNumber = parameterSet.get<unsigned int>("arbitraryElementNumber", 0);
-//
-// FieldVector<int,3> row;
+//
+// FieldVector<int,3> row;
// row = arbitraryComponentwiseIndices(Basis_CE,arbitraryElementNumber,arbitraryLeafIndex);
// printvector(std::cout, row, "row" , "--");
-//
+//
@@ -1704,12 +1704,12 @@ int main(int argc, char *argv[])
VectorCT x_1 = load_alpha1;
VectorCT x_2 = load_alpha2;
VectorCT x_3 = load_alpha3;
-
+
auto load_alpha1BS = load_alpha1;
-
+
// printvector(std::cout, load_alpha1, "load_alpha1 before SOLVER", "--" );
// printvector(std::cout, load_alpha2, "load_alpha2 before SOLVER", "--" );
-
+
if (Solvertype == 1) // CG - SOLVER
{
std::cout << "------------ CG - Solver ------------" << std::endl;
@@ -1735,10 +1735,10 @@ int main(int argc, char *argv[])
solver.apply(x_3, load_alpha3, statistics);
}
////////////////////////////////////////////////////////////////////////////////////
-
+
else if (Solvertype ==2) // GMRES - SOLVER
{
-
+
std::cout << "------------ GMRES - Solver ------------" << std::endl;
// Turn the matrix into a linear operator
MatrixAdapter<MatrixCT,VectorCT,VectorCT> stiffnessOperator(stiffnessMatrix_CE);
@@ -1783,12 +1783,12 @@ int main(int argc, char *argv[])
// printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
// printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
// printvector(std::cout, load_alpha2, "load_alpha2 AFTER SOLVER", "--" );
-
+
////////////////////////////////////////////////////////////////////////////////////
// Extract phi_alpha & M_alpha coefficients
////////////////////////////////////////////////////////////////////////////////////
-
+
VectorCT phi_1, phi_2, phi_3;
phi_1.resize(Basis_CE.size());
phi_1 = 0;
@@ -1804,16 +1804,16 @@ int main(int argc, char *argv[])
phi_2[i] = x_2[i];
phi_3[i] = x_3[i];
}
-
+
FieldVector<double,3> m_1, m_2, m_3;
-
+
for(size_t i=0; i<3; i++)
{
m_1[i] = x_1[phiOffset+i];
m_2[i] = x_2[phiOffset+i];
m_3[i] = x_3[phiOffset+i];
}
-
+
// assemble M_alpha's (actually iota(M_alpha) )
MatrixRT M_1(0), M_2(0), M_3(0);
@@ -1838,7 +1838,7 @@ int main(int argc, char *argv[])
////////////////////////////////////////////////////////////////////////////
// Substract Integral-mean
- //////////////////////////////////////////////////////////////////////////// // ERROR
+ //////////////////////////////////////////////////////////////////////////// // ERROR
using SolutionRange = FieldVector<double,dim>;
if(write_IntegralMean)
@@ -1856,7 +1856,7 @@ int main(int argc, char *argv[])
subtractIntegralMean(Basis_CE, phi_1);
subtractIntegralMean(Basis_CE, phi_2);
subtractIntegralMean(Basis_CE, phi_3);
- subtractIntegralMean(Basis_CE, x_1);
+ subtractIntegralMean(Basis_CE, x_1);
subtractIntegralMean(Basis_CE, x_2);
subtractIntegralMean(Basis_CE, x_3);
//////////////////////////////////////////
@@ -1874,29 +1874,29 @@ int main(int argc, char *argv[])
/////////////////////////////////////////////////////////
// Write Solution (Corrector Coefficients) in Logs
- /////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////
log << "\nSolution of Corrector problems:\n";
if(write_corrector_phi1)
{
- log << "\n Corrector_phi1:\n";
+ log << "\n Corrector_phi1:\n";
log << x_1 << std::endl;
}
if(write_corrector_phi2)
{
log << "-----------------------------------------------------";
- log << "\n Corrector_phi2:\n";
+ log << "\n Corrector_phi2:\n";
log << x_2 << std::endl;
}
if(write_corrector_phi3)
{
log << "-----------------------------------------------------";
- log << "\n Corrector_phi3:\n";
+ log << "\n Corrector_phi3:\n";
log << x_3 << std::endl;
}
////////////////////////////////////////////////////////////////////////////
// Make a discrete function from the FE basis and the coefficient vector
- //////////////////////////////////////////////////////////////////////////// // ERROR
+ //////////////////////////////////////////////////////////////////////////// // ERROR
using SolutionRange = FieldVector<double,dim>;
auto correctorFunction_1 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
@@ -1909,14 +1909,14 @@ int main(int argc, char *argv[])
auto correctorFunction_3 = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(
Basis_CE,
phi_3);
-
+
/////////////////////////////////////////////////////////
- // Create Containers for Basis-Matrices, Correctors and Coefficients
+ // Create Containers for Basis-Matrices, Correctors and Coefficients
/////////////////////////////////////////////////////////
const std::array<Func2Tensor, 3> x3MatrixBasis = {x3G_1, x3G_2, x3G_3};
- const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
+ const std::array<VectorCT, 3> coeffContainer = {x_1, x_2, x_3};
const std::array<VectorCT, 3> loadContainer = {load_alpha1, load_alpha2, load_alpha3};
const std::array<MatrixRT*, 3 > mContainer = {&M_1 , &M_2, & M_3};
@@ -1929,20 +1929,20 @@ int main(int argc, char *argv[])
VectorCT tmp1,tmp2;
FieldVector<double,3> B_hat ;
- // compute Q
+ // compute Q
for(size_t a = 0; a < 3; a++)
for(size_t b=0; b < 3; b++)
{
assembleCellLoad(Basis_CE, muLocal, lambdaLocal, gamma, tmp1 ,x3MatrixBasis[b]); // <L i(M_alpha) + sym(grad phi_alpha), i(x3G_beta) >
double GGterm = 0.0;
- GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
+ GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , x3MatrixBasis[b] ); // <L i(x3G_alpha) , i(x3G_beta) >
Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric... check positiv definitness?
-// Q[a][b] = (coeffContainer[a]*loadContainer[b]) + GGterm; // TEST
+// Q[a][b] = (coeffContainer[a]*loadContainer[b]) + GGterm; // TEST
std::cout << "GGTerm:" << GGterm << std::endl;
std::cout << "coeff*tmp: " << coeffContainer[a]*tmp1 << std::endl;
-// std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
+// std::cout << "Test : " << coeffContainer[a]*loadContainer[b] << std::endl; //same...
}
printmatrix(std::cout, Q, "Matrix Q", "--");
log << "Computed Matrix Q: " << std::endl;
@@ -1956,6 +1956,7 @@ int main(int argc, char *argv[])
auto GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3MatrixBasis[a] , B_Term); // <L i(x3G_alpha) , B>
B_hat[a] = (coeffContainer[a]*tmp2) + GGterm;
+ std::cout << "check this Contribution: " << (coeffContainer[a]*tmp2) << std::endl; //see orthotropic.tex
std::cout <<"B_hat[i]: " << B_hat[a] << std::endl;
}
log << "Computed prestrain B_hat: " << std::endl;
@@ -1978,7 +1979,7 @@ int main(int argc, char *argv[])
-
+
auto q1_c = Q[0][0];
auto q2_c = Q[1][1];
auto q3_c = Q[2][2];
@@ -1997,11 +1998,11 @@ int main(int argc, char *argv[])
log << "computed b3_hat: " << B_hat[2] << std::endl;
-
+
//////////////////////////////////////////////////////////////
// Define Analytic Solutions
//////////////////////////////////////////////////////////////
-
+
// double b1 = (mu1*p1/4)*(beta/(theta+(1-theta)*beta))*(1-theta*(1+alpha));
// double b2 = (mu1*p1/8)*(1-theta*(1+beta*alpha));
@@ -2026,42 +2027,45 @@ int main(int argc, char *argv[])
log << "q2 : " << q2 << std::endl;
log << "q3 should be between q1 and q2" << std::endl;
-
+
auto symPhi_1_analytic = [mu1, mu2, theta, muTerm] (const Domain& x) {
return MatrixRT{{ ((mu1*mu2)/((theta*mu1 +(1-theta)*mu2)*muTerm(x)) - 1)*x[2] , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
-
+
auto zeroFunction = [] (const Domain& x) {
return MatrixRT{{0.0 , 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
-
-
+
+
if(write_L2Error)
{
- auto L2error = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- std::cout << "L2-Error: " << L2error << std::endl;
-
+ auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ std::cout << "L2-SymError: " << L2SymError << std::endl;
+
// auto L2errorTest = computeL2ErrorTest(Basis_CE,phi_1,gamma,symPhi_1_analytic);
// std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl;
-
+
auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl; // TODO Not Needed
-
+
VectorCT zeroVec;
zeroVec.resize(Basis_CE.size());
zeroVec = 0;
auto L2Norm_analytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
std::cout << "L2-Norm(analytic): " << L2Norm_analytic << std::endl;
-
- log << "L2-Error (symmetric Gradient phi_1):" << L2error << std::endl;
+
+ log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
+
+ //TEST
+ std::cout << "L2ErrorOld:" << computeL2ErrorOld(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
}
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////
-
+
VTKWriter<GridView> vtkWriter(gridView_CE);
-
+
vtkWriter.addVertexData(
correctorFunction_1,
VTK::FieldInfo("corrector phi_1", VTK::FieldInfo::Type::vector, dim));
@@ -2073,6 +2077,6 @@ int main(int argc, char *argv[])
VTK::FieldInfo("corrector phi_3", VTK::FieldInfo::Type::vector, dim));
vtkWriter.write("CellProblem-result");
std::cout << "wrote data to file: CellProblem-result" << std::endl; // better with string for output name..
-
+
log.close();
}