diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 0fb19495532575032af50a86c03fa3881117d8a8..dd89c51063406e8a42e57220a70e988abe233285 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -14,15 +14,26 @@ cellDomain = 1
# Grid parameters
#############################################
-#nElements_Cell = 20 20 20 # number elements in each direction (y1 y2 x3)
+levels = 2
+
+#
+#Elements_Cell = 20 20 20 # number elements in each direction (y1 y2 x3)
+#nElements_Cell = 30 30 30
+#nElements_Cell = 30 30 30
+#nElements_Cell = 50 50 50
#nElements_Cell = 100 100 2
-nElements_Cell = 10 10 10
+#nElements_Cell = 100 100 100 // does not work
+#nElements_Cell = 10 10 10
+
#nElements_Cell = 4 4 4
+#nElements_Cell = 8 8 8
+#nElements_Cell = 16 16 16
+nElements_Cell = 32 32 32
width = 1.0
-gamma = 50.0
+gamma = 1.0
#############################################
# Material parameters
@@ -83,9 +94,13 @@ Func2Tensor B2_ = [this] (const Domain& x) { // ISOTROPIC PRESSURE
#############################################
set_IntegralZero = true
+#set_IntegralZero = false
+
+#arbitraryLocalIndex = 7
+#arbitraryElementNumber = 3
-arbitraryLocalIndex = 7
-arbitraryElementNumber = 3
+arbitraryLocalIndex = 0
+arbitraryElementNumber = 0
#############################################
@@ -95,9 +110,12 @@ arbitraryElementNumber = 3
Solvertype = 1
-write_corrector_phi1 = false
-write_corrector_phi2 = false
-write_corrector_phi3 = false
+#write_corrector_phi1 = false
+#write_corrector_phi2 = false
+#write_corrector_phi3 = false
+#write_corrector_phi1 = true
+#write_corrector_phi2 = true
+#write_corrector_phi3 = true
write_L2Error = true
write_IntegralMean = true
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 31a59ba66eb6e414ed4ea52cb6dde193d06d5d04..58d7d1deac04a0a835720ab87d61ed9ebe7e64dd 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -820,7 +820,7 @@ auto energy(const Basis& basis,
energyDensity += stressU[ii][jj] * strain2[ii][jj];
// elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
- elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+ elementEnergy += energyDensity * quadPoint.weight() * integrationElement; //TODO integratonElement needed here?
// elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
@@ -1009,94 +1009,6 @@ 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,
- Vector& coeffVector,
- const double gamma,
- Domain& x // evaluation Point in reference coordinates
- )
-{
-
-
- constexpr int dim = 3;
- constexpr int nCompo = 3;
-
- auto localView = basis.localView();
-
-
- using GridView = typename Basis::GridView;
-// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
-
-
-
-
-
- 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,
- referenceGradients);
-
- // Compute the shape function gradients on the grid element
- std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
- // std::vector< VectorRT> gradients(referenceGradients.size());
-
- for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-
- 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 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++)
- for (int jj=0; jj<nCompo; jj++)
- {
- strainU[ii][jj] = coeffVector[globalIdx] * 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // symmetric gradient
- }
-
-
- }
- }
-
-}*/
-
-
-
template<class Basis, class Vector, class MatrixFunction>
double computeL2SymErrorNew(const Basis& basis,
Vector& coeffVector,
@@ -1150,11 +1062,11 @@ double computeL2SymErrorNew(const Basis& basis,
MatrixRT tmp(0);
-
+ MatrixRT tmpNew(0);
double sum = 0.0;
- for (size_t i=0; i < nSf; i++)
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); // hier i:leafIdx
@@ -1166,18 +1078,22 @@ double computeL2SymErrorNew(const Basis& basis,
defGradientU[k][1] = coeffVector[globalIdx1]*gradients[i][1]; //X2
defGradientU[k][2] = coeffVector[globalIdx1]*(1.0/gamma)*gradients[i][2]; //X3
+// printvector(std::cout, gradients[i], "gradients[i]", "--");
+// std::cout << "coeffVector[globalIdx1]" << coeffVector[globalIdx1] << std::endl;
+
// 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]);
- }
+// 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]);
+// }
-// printmatrix(std::cout, strainU, "strainU", "--");
-// printmatrix(std::cout, tmp, "tmp", "--");
- tmp += strainU;
+// printmatrix(std::cout, strainU, "strainU", "--");
// printmatrix(std::cout, tmp, "tmp", "--");
+// tmp += strainU;
+ tmpNew += defGradientU;
+// printmatrix(std::cout, tmp, "tmp", "--");
// Local energy density: stress times strain
@@ -1191,14 +1107,33 @@ double computeL2SymErrorNew(const Basis& basis,
// }
}
- tmp = tmp - matrixField(quadPos);
+// printmatrix(std::cout, tmpNew, "tmpNew", "--");
for (int ii=0; ii<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
- sum += std::pow(tmp[ii][jj],2);
+ {
+ tmp[ii][jj] = 0.5 * (tmpNew[ii][jj] + tmpNew[jj][ii]);
+ }
- error += sum * quadPoint.weight() * integrationElement;
+// printmatrix(std::cout, matrixField(quadPos), "matrixField(quadPos)", "--");
+// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat ganz andere Struktur als analytic?
+
+
+// tmp = tmp - matrixField(quadPos);
+// printmatrix(std::cout, tmp - matrixField(quadPos), "Difference", "--");
+
+
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ {
+ sum += std::pow(tmp[ii][jj]-matrixField(quadPos)[ii][jj],2);
+// sum += std::pow(tmp[ii][jj],2);
+ }
+// std::cout << "sum:" << sum << std::endl;
+
+ error += sum * quadPoint.weight() * integrationElement;
+// std::cout << "error:" << error << std::endl;
}
}
@@ -2187,8 +2122,18 @@ int main(int argc, char *argv[])
FieldVector<double,dim> lower({0.0, 0.0, -1.0/2.0});
FieldVector<double,dim> upper({1.0, 1.0, 1.0/2.0});
}
+
+
+
std::array<int, dim> nElements = parameterSet.get<std::array<int,dim>>("nElements_Cell", {10, 10, 10});
+
+
+ ///// LEVEL - APPROACH // TODO
+/* int levels = parameterSet.get<int>("levels", 1); //besser : numLevels
+ std::array<int, dim> nElements = { std::pow(2,levels) , std::pow(2,levels) , std::pow(2,levels) }; */
+
+ std::cout << "Number of Elements in each direction: " << nElements << std::endl;
log << "Number of Elements in each direction: " << nElements << std::endl;
using CellGridType = YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> >;
@@ -2220,10 +2165,10 @@ int main(int argc, char *argv[])
// Create Lambda-Functions for material Parameters depending on microstructure
///////////////////////////////////
auto muTerm = [mu1, mu2, theta] (const Domain& z) {
-// if (abs(z[0]) >= (theta/2.0)) //TODO check ..nullset/boundary
-// return mu1;
- if (abs(z[0]) > (theta/2.0)) //TEST include boundary (indicatorFunction)
+ if (abs(z[0]) >= (theta/2.0)) //TODO check ..nullset/boundary
return mu1;
+// if (abs(z[0]) > (theta/2.0)) //TEST include boundary (indicatorFunction)
+// return mu1;
// if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
// return mu2;
else
@@ -2780,7 +2725,7 @@ int main(int argc, char *argv[])
std::cout << " ----- TEST ----" << std::endl;
- std::cout << "computeL2SymErrorNew: " << computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << "computeL2SymErrorNew: " << computeL2SymErrorNew(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
std::cout << " -----------------" << std::endl;
auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
@@ -2793,7 +2738,7 @@ int main(int argc, char *argv[])
std::cout << "L2-Norm(Symphi_1) w zerofunction: " << L2Norm_Symphi << std::endl; // TODO Not Needed
- std::cout << "L2 -Norm(SymPhi_1): " << computeL2SymNormCoeff(Basis_CE,phi_1,gamma) << std::endl; // does not make a difference
+ std::cout << "L2-Norm(SymPhi_1): " << computeL2SymNormCoeff(Basis_CE,phi_1,gamma) << std::endl; // does not make a difference
VectorCT zeroVec;
zeroVec.resize(Basis_CE.size());