diff --git a/inputs/cellsolver.parset b/inputs/cellsolver.parset
index 1e700a15b1db5e2ea7a4aa0c3f4c180d56587359..19b4e6848d73f8d7a11083a2c3d23dc8165cf58e 100644
--- a/inputs/cellsolver.parset
+++ b/inputs/cellsolver.parset
@@ -3,6 +3,14 @@
#path for logfile
outputPath = "../../outputs/output.txt"
+
+#############################################
+# Debug Output
+#############################################
+
+#print_debug = true #default = false
+
+
#############################################
# Cell Domain
#############################################
@@ -20,9 +28,9 @@ cellDomain = 1
## {start,finish} computes on all grid from 2^(start) to 2^finish refinement
########################################################################
-numLevels = 1 4 # computes all levels from first to second entry
+numLevels = 5 5 # computes all levels from first to second entry
#numLevels = 3 3 # computes all levels from first to second entry
-
+#numLevels = 1 6
@@ -46,8 +54,8 @@ numLevels = 1 4 # computes all levels from first to second entry
width = 1.0
-#gamma = 50.0
-gamma = 1.0
+gamma = 50.0
+#gamma = 1.0
#############################################
# Material parameters
@@ -63,8 +71,8 @@ lambda1 = 0.0
# Prestrain parameters
#############################################
-#theta = 0.3 # volume fraction
-theta = 0.25 # volume fraction
+#theta = 0.3 # volume fraction #default = 1.0/4.0
+#theta = 0.25 # volume fraction
#theta = 0.75 # volume fraction
prestrainType = 2
@@ -134,7 +142,7 @@ Solvertype = 1
#write_corrector_phi3 = true
write_L2Error = true
-write_IntegralMean = true
+#write_IntegralMean = true
#############################################
# Define Analytic Solutions
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index c3b7f72cd4ecd6e2baca3b52d51850f5585696fb..08911b48fe0b24019942ef0e05d8bb21d5e13587 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -868,7 +868,7 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
template<class Basis, class Vector, class MatrixFunction>
-double computeL2SymErrorNew(const Basis& basis,
+double computeL2SymError(const Basis& basis,
Vector& coeffVector,
const double gamma,
const MatrixFunction& matrixFieldFunc
@@ -879,10 +879,8 @@ double computeL2SymErrorNew(const Basis& basis,
constexpr int dimWorld = 3;
auto localView = basis.localView();
-
auto matrixFieldGVF = Dune::Functions::makeGridViewFunction(matrixFieldFunc, basis.gridView());
auto matrixField = localFunction(matrixFieldGVF);
-
using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
for (const auto& element : elements(basis.gridView()))
@@ -894,21 +892,20 @@ double computeL2SymErrorNew(const Basis& basis,
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
const auto nSf = localFiniteElement.localBasis().size();
- int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1 ); //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();
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPos);
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
+ const auto integrationElement = geometry.integrationElement(quadPos);
- std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPos,
- referenceGradients);
+ std::vector< FieldMatrix<double, 1, dim>> referenceGradients;
+ localFiniteElement.localBasis().evaluateJacobian(quadPos,referenceGradients);
// Compute the shape function gradients on the grid element
- std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
+ std::vector<FieldVector<double,dim>> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
@@ -919,7 +916,6 @@ double computeL2SymErrorNew(const Basis& basis,
for (size_t k=0; k < dimWorld; k++)
for (size_t i=0; i < nSf; i++)
{
-
size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
size_t globalIdx1 = localView.index(localIdx1);
@@ -930,23 +926,22 @@ double computeL2SymErrorNew(const Basis& basis,
defGradientU[k][2] = coeffVector[globalIdx1]*(1.0/gamma)*gradients[i][2]; //X3
// symmetric Gradient (Elastic Strains)
- MatrixRT strainU(0);
+// MatrixRT strainU(0);
for (int ii=0; ii<dimWorld; ii++)
for (int jj=0; jj<dimWorld; jj++)
{
- strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
+// strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]); // not needed (old version)
+ tmp[ii][jj] += 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
+// tmp += strainU;
// printmatrix(std::cout, strainU, "strainU", "--");
// printmatrix(std::cout, tmp, "tmp", "--");
- tmp += strainU;
+//
}
-
// printmatrix(std::cout, matrixField(quadPos), "matrixField(quadPos)", "--");
-// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat ganz andere Struktur als analytic?
-
+// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat andere Struktur als analytic?
// tmp = tmp - matrixField(quadPos);
// printmatrix(std::cout, tmp - matrixField(quadPos), "Difference", "--");
-
for (int ii=0; ii<dimWorld; ii++)
for (int jj=0; jj<dimWorld; jj++)
{
@@ -962,589 +957,9 @@ double computeL2SymErrorNew(const Basis& basis,
-
-
-
-
-template<class Basis, class Vector, class MatrixFunction>
-double computeL2SymError(const Basis& basis,
- Vector& coeffVector,
- const double gamma,
- const MatrixFunction& matrixFieldFunc
- )
-{
- // This Version computes the SCALAR PRODUCT of the difference ..
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
-
- auto localView = basis.localView();
-
- auto matrixFieldGVF = Dune::Functions::makeGridViewFunction(matrixFieldFunc, basis.gridView());
- auto matrixField = localFunction(matrixFieldGVF);
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
- 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();
- 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);
-
-
- for (const auto& quadPoint : quad)
- {
- const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
- 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), defGradientV(0);
-
-
- for (size_t k=0; k < dimWorld; k++)
- for (size_t i=0; i < nSf; i++)
- {
- for (size_t l=0; l < dimWorld; l++)
- for (size_t j=0; j < nSf; j++ )
- {
-
- 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 globalIdx2 = localView.index(localIdx2);
-
- MatrixRT defGradientU(0);
- // (scaled) Deformation gradient of the ansatz basis function
- 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
-
- MatrixRT defGradientV(0);
- defGradientV[l][0] = gradients[j][0]; // Y
- defGradientV[l][1] = gradients[j][1]; //X2
- defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
-
- // symmetric Gradient (Elastic Strains)
- MatrixRT strainU(0), strainV(0);
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- {
- 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<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
-
- error += tmp * quadPoint.weight() * integrationElement;
- }
- }
-
-
- for (size_t k=0; k < dimWorld; k++)
- for (size_t i=0; i < nSf; i++)
- {
- size_t localIdx1 = localView.tree().child(k).localIndex(i);
- size_t globalIdx1 = localView.index(localIdx1);
-
- // (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(0);
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- {
- strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
- }
-
- double tmp = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; 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<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- tmp += matrixField(quadPos)[ii][jj] * matrixField(quadPos)[ii][jj];
-
- error += tmp * quadPoint.weight() * integrationElement;
-
- }
- }
- return sqrt(error);
-}
-
-
-
- ////////////////////////////////////////////////////////////// L2-NORM of symmetric Gradient /////////////////////////////////////////////////////////////////
-
-
+////////////////////////////////////////////////////////////// L2-NORM /////////////////////////////////////////////////////////////////
template<class Basis, class Vector>
-double computeL2SymNormCoeff(const Basis& basis,
- Vector& coeffVector,
- const double gamma
- )
-{
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
-
- auto localView = basis.localView();
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
- for (const auto& element : elements(basis.gridView()))
- {
- localView.bind(element);
- auto geometry = element.geometry();
-
- const auto& localFiniteElement = localView.tree().child(0).finiteElement();
- 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);
-
- for (const auto& quadPoint : quad)
- {
-
- const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
- referenceGradients);
-
- // Compute the shape function gradients on the grid element
- std::vector<FieldVector<double,dim> > gradients(referenceGradients.size());
-
- for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
-
-
- MatrixRT tmp(0);
- double sum = 0.0;
-
- for (size_t i=0; i < nSf; i++)
- for (size_t k=0; k < dimWorld; k++)
- {
-
- size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
- size_t globalIdx1 = localView.index(localIdx1);
-
- // (scaled) Deformation gradient of the ansatz basis function
- MatrixRT defGradientU(0);
- defGradientU[k][0] = coeffVector[globalIdx1]*gradients[i][0]; // Y //hier i:leaf oder localIdx?
- defGradientU[k][1] = coeffVector[globalIdx1]*gradients[i][1]; //X2
- defGradientU[k][2] = coeffVector[globalIdx1]*(1.0/gamma)*gradients[i][2]; //X3
-
- // symmetric Gradient (Elastic Strains)
- MatrixRT strainU(0);
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; 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, tmp, "tmp", "--");
-
-
- // Local energy density: stress times strain
-// double tmp = 0.0;
-// for (int ii=0; ii<dimWorld; ii++)
-// for (int jj=0; jj<dimWorld; jj++)
-// tmp[ii][jj] += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
-//
-//
-// error += tmp * quadPoint.weight() * integrationElement;
-// }
- }
-
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- sum += std::pow(tmp[ii][jj],2);
-
- error += sum * quadPoint.weight() * integrationElement;
-
- }
- }
-
- return sqrt(error);
-}
-
-
-
-
-
-
-
-
-
-// TEST
-template<class Basis, class Vector>
-double computeL2SymNormCoeffV2(const Basis& basis,
- Vector& coeffVector,
- const double gamma
- )
-{
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
-
- auto localView = basis.localView();
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
- 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();
-// 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);
-
- for (const auto& quadPoint : quad)
- {
-
- const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
- 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), defGradientV(0);
-
-
- for (size_t k=0; k < dimWorld; k++)
- for (size_t i=0; i < nSf; i++)
- {
- for (size_t l=0; l < dimWorld; l++)
- for (size_t j=0; j < nSf; j++ )
- {
-
- 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 globalIdx2 = localView.index(localIdx2);
-
- // (scaled) Deformation gradient of the ansatz basis function
- MatrixRT defGradientU(0);
- 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
-
- MatrixRT defGradientV(0);
- defGradientV[l][0] = gradients[j][0]; // Y
- defGradientV[l][1] = gradients[j][1]; //X2
- defGradientV[l][2] = (1.0/gamma)*gradients[j][2]; //X3
-
- // symmetric Gradient (Elastic Strains)
- MatrixRT strainU(0), strainV(0);
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- {
- 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<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- tmp += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
-
-
- error += tmp * quadPoint.weight() * integrationElement;
- }
- }
-
-
- }
- }
-
- return sqrt(error);
-}
-
-
-
-
-
-
-
-
- ////////////////////////////////////////////////////////////// L2-NORM /////////////////////////////////////////////////////////////////
-// TEST
-template<class Basis, class Vector>
-double computeL2NormCoeff(const Basis& basis,
- Vector& coeffVector
- )
-{
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
- auto localView = basis.localView();
-
- 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();
-// 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);
-
-
- for (const auto& quadPoint : quad)
- {
- const auto& quadPos = quadPoint.position();
-// const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldVector<double, 1>> functionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPoint.position(),
- functionValues);
-
-// for (auto&& value : functionValues)
-// std::cout << value << std::endl;
-
- /*
- using LocalFEType = LagrangeSimplexLocalFiniteElement<double,double,dim,2>;
- LocalFEType localFE;
- // Get shape function values
- using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
- std::vector<ValueType> values;*/
-
- double tmp = 0.0;
- for (size_t k=0; k < dimWorld; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
- {
- double comp = 0.0;
- for (size_t i=0; i < nSf; i++)
- {
- // for (size_t l=0; l < dimWorld; l++)
- // for (size_t j=0; j < nSf; j++ )
- // {
- size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
- size_t globalIdx1 = localView.index(localIdx1);
- // size_t localIdx2 = localView.tree().child(k).localIndex(j);
- // size_t localIdx2 = localView.tree().child(l).localIndex(j);
- // size_t globalIdx2 = localView.index(localIdx2);
-
-
- comp += coeffVector[globalIdx1]*functionValues[i];
-
- // tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
- // tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
- // std::cout << "tmp:" << tmp << std::endl;
-
-// error += tmp * quadPoint.weight() * integrationElement;
- }
- tmp += std::pow(comp,2);
- }
- error += tmp * quadPoint.weight() * integrationElement;
- }
- }
- return sqrt(error);
-}
-
-
-// TEST
-template<class Basis, class Vector>
-double computeL2NormCoeffV2(const Basis& basis,
- Vector& coeffVector
- )
-{
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
- auto localView = basis.localView();
-
-// using GridView = typename Basis::GridView;
-// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-// using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
- 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();
-// 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);
-
-
- for (const auto& quadPoint : quad)
- {
- const auto& quadPos = quadPoint.position();
-// const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldVector<double, 1>> functionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPoint.position(),
- functionValues);
-
-// for (auto&& value : functionValues)
-// std::cout << value << std::endl;
-
- /*
- using LocalFEType = LagrangeSimplexLocalFiniteElement<double,double,dim,2>;
- LocalFEType localFE;
- // Get shape function values
- using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
- std::vector<ValueType> values;*/
-
- for (size_t k=0; k < dimWorld; k++) // ANALOG STOKES Beitrag nur wenn k=l!!!
- for (size_t i=0; i < nSf; i++)
- {
-// for (size_t l=0; l < dimWorld; l++)
- for (size_t j=0; j < nSf; j++ )
- {
- size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
- size_t globalIdx1 = localView.index(localIdx1);
- size_t localIdx2 = localView.tree().child(k).localIndex(j);
-// size_t localIdx2 = localView.tree().child(l).localIndex(j);
- size_t globalIdx2 = localView.index(localIdx2);
-
- double tmp = 0.0;
- tmp += coeffVector[globalIdx1]*functionValues[i]*coeffVector[globalIdx2]*functionValues[j];
-
-// tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
-// tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
-// std::cout << "tmp:" << tmp << std::endl;
-
- error += tmp * quadPoint.weight() * integrationElement;
- }
- }
- }
-
- }
- return sqrt(error);
-}
-
-// // TEST
-template<class Basis, class Vector>
-double computeL2NormCoeffV3(const Basis& basis, // Falsch: betrachtet keine gemischten TERME!
- Vector& coeffVector
- )
-{
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
- auto localView = basis.localView();
-
-// using GridView = typename Basis::GridView;
-// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
-// using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
- 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();
-// 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);
-
-
- for (const auto& quadPoint : quad)
- {
-
- const auto& quadPos = quadPoint.position();
-// const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
-
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldVector<double, 1>> functionValues;
- localFiniteElement.localBasis().evaluateFunction(quadPoint.position(),
- functionValues);
-
-// for (auto&& value : functionValues)
-// std::cout << value << std::endl;
-
- /*
- using LocalFEType = LagrangeSimplexLocalFiniteElement<double,double,dim,2>;
- LocalFEType localFE;
- // Get shape function values
- using ValueType = LocalFEType::Traits::LocalBasisType::Traits::RangeType;
- std::vector<ValueType> values;*/
-
- for (size_t k=0; k < dimWorld; k++)
- for (size_t i=0; i < nSf; i++)
- {
- size_t localIdx1 = localView.tree().child(k).localIndex(i); // hier i:leafIdx
- size_t globalIdx1 = localView.index(localIdx1);
-
- double tmp = 0.0;
-
- tmp = std::pow(coeffVector[globalIdx1]*functionValues[i],2); // same value for each component k...
-// std::cout << "tmp:" << tmp << std::endl;
-
- error += tmp * quadPoint.weight() * integrationElement;
- }
- }
- }
- return sqrt(error);
-}
-
-
-// TEST
-template<class Basis, class Vector>
-double computeL2NormFunc(const Basis& basis,
+double computeL2Norm(const Basis& basis,
Vector& coeffVector
)
{
@@ -1553,17 +968,9 @@ double computeL2NormFunc(const Basis& basis,
constexpr int dim = 3;
constexpr int dimWorld = 3;
auto localView = basis.localView();
-
-
- using SolutionRange = FieldVector<double,dim>;
- auto GVFunc = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(basis,coeffVector);
+ auto GVFunc = Functions::makeDiscreteGlobalBasisFunction<FieldVector<double,dim>>(basis,coeffVector);
auto localfun = localFunction(GVFunc);
-
-// FieldVector<double,3> r = {0.0, 0.0, 0.0};
-// double r = 0.0;
-
- // Compute Area integral & integral of FE-function
for(const auto& element : elements(basis.gridView()))
{
localView.bind(element);
@@ -1575,117 +982,18 @@ double computeL2NormFunc(const Basis& basis,
for(const auto& quadPoint : quad)
{
- const double integrationElement = element.geometry().integrationElement(quadPoint.position());
-
const auto& quadPos = quadPoint.position();
+ const double integrationElement = element.geometry().integrationElement(quadPos);
double tmp = localfun(quadPos) * localfun(quadPos);
error += tmp * quadPoint.weight() * integrationElement;
-
}
}
-
- // std::cout << "Domain-Area: " << area << std::endl;
return sqrt(error);
}
-
-
-
-template<class Basis, class Vector, class MatrixFunction>
-double computeL2ErrorOld(const Basis& basis,
- Vector& coeffVector,
- const double gamma,
- const MatrixFunction& matrixFieldFunc
- )
-{
-
- auto error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
-
- auto localView = basis.localView();
-
- auto matrixFieldGVF = Dune::Functions::makeGridViewFunction(matrixFieldFunc, basis.gridView());
- auto matrixField = localFunction(matrixFieldGVF);
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
-
- 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);
- const auto& quad = QuadratureRules<double,dim>::rule(element.type(), orderQR);
-
-
- for (const auto& quadPoint : quad)
- {
- const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
- 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 < dimWorld; k++)
- for (size_t i=0; i < nSf; 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] = 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<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; 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
- }
-
- // Local energy density: stress times strain
- double tmp = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; 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);
-}
-
+// OPTIONAL H1-Seminorm ...
/*
template<class VectorCT>
@@ -1756,104 +1064,6 @@ for (const auto& e : elements(basis_.gridView()))
-template<class Basis, class Vector, class MatrixFunction>
-double computeL2SymErrorSPTest(const Basis& basis,
- Vector& coeffVector,
- const double gamma,
- const MatrixFunction& matrixFieldFunc
- )
-{
- // This Version computes the SCALAR PRODUCT of the difference ..
- double error = 0.0;
- constexpr int dim = 3;
- constexpr int dimWorld = 3;
-
- auto localView = basis.localView();
-
- auto matrixFieldGVF = Dune::Functions::makeGridViewFunction(matrixFieldFunc, basis.gridView());
- auto matrixField = localFunction(matrixFieldGVF);
-
- using GridView = typename Basis::GridView;
- using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
- using MatrixRT = FieldMatrix< double, dimWorld, dimWorld>;
-
- Vector zeroVec;
- zeroVec.resize(basis.size());
- zeroVec = 0;
-
- double Term1 = std::pow(computeL2SymNormCoeff(basis,coeffVector,gamma),2);
-// std::cout << "Term1:" << Term1 << std::endl;
- double Term2 = std::pow(computeL2SymError(basis,zeroVec ,gamma,matrixFieldFunc) ,2);
-// std::cout << "Term2:" << Term2 << std::endl;
-
-
- 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();
-// 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);
-
-
-
-
- for (const auto& quadPoint : quad)
- {
- const auto& quadPos = quadPoint.position();
- const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
- const auto integrationElement = geometry.integrationElement(quadPoint.position());
-
- std::vector< FieldMatrix< double, 1, dim> > referenceGradients;
- localFiniteElement.localBasis().evaluateJacobian(quadPoint.position(),
- 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]); // TRANSFORMED
-
-
- for (size_t k=0; k < dimWorld; k++)
- for (size_t i=0; i < nSf; i++)
- {
- size_t localIdx1 = localView.tree().child(k).localIndex(i);
- size_t globalIdx1 = localView.index(localIdx1);
-
- // (scaled) Deformation gradient of the ansatz basis function
- MatrixRT defGradientU(0); //TEST (doesnt change anything..)
- 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<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- {
- strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
- }
-
- double tmp = 0.0;
- for (int ii=0; ii<dimWorld; ii++)
- for (int jj=0; jj<dimWorld; jj++)
- tmp += (-2.0)*coeffVector[globalIdx1]*strainU[ii][jj] * matrixField(quadPos)[ii][jj];
-
- error += tmp * quadPoint.weight() * integrationElement;
- }
-
- }
- }
-
-
- return sqrt(error+Term1+Term2);
-}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -2549,49 +1759,34 @@ int main(int argc, char *argv[])
if(write_L2Error)
{
- std::cout << " ----- TEST L2ErrorSym ----" << std::endl;
- std::cout << "computeL2SymErrorNew: " << computeL2SymErrorNew(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << " ----- L2ErrorSym ----" << std::endl;
auto L2SymError = computeL2SymError(Basis_CE,phi_1,gamma,symPhi_1_analytic);
- std::cout << "L2-SymError: " << L2SymError << std::endl;
- std::cout << "computeL2SymErrorSPTest: " << computeL2SymErrorSPTest(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << "L2SymError: " << L2SymError << std::endl;
std::cout << " -----------------" << std::endl;
-
-
- std::cout << " ----- TEST L2NORMSym ----" << std::endl;
- auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction);
- 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 << " ----- L2NORMSym ----" << std::endl;
+ auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction);
+ std::cout << "L2-Norm(Symphi_1): " << L2Norm_Symphi << std::endl;
VectorCT zeroVec;
zeroVec.resize(Basis_CE.size());
zeroVec = 0;
- auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
- std::cout << "L2-Norm(Symanalytic): " << L2Norm_SymAnalytic << std::endl;
-
+ auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma, symPhi_1_analytic);
+ std::cout << "L2-Norm(SymAnalytic): " << L2Norm_SymAnalytic << std::endl;
std::cout << " -----------------" << std::endl;
-
- std::cout << " ----- TEST L2NORM ----" << std::endl;
- auto L2Norm = computeL2NormCoeff(Basis_CE,phi_1);
- std::cout << "L2Norm(phi_1) - coeff: " << L2Norm << std::endl;
-// std::cout << "L2Norm(phi_1) - coeff: " << computeL2NormCoeff(Basis_CE,phi_1) << std::endl;
- std::cout << "L2Norm(phi_1) - GVfunc: " << computeL2NormFunc(Basis_CE,phi_1) << std::endl;
- std::cout << "L2Norm(phi_1) - coeffV2: " << computeL2NormCoeffV2(Basis_CE,phi_1) << std::endl;
- std::cout << "L2Norm(phi_1) - coeffV3: " << computeL2NormCoeffV3(Basis_CE,phi_1) << std::endl;
- std::cout << "L2ErrorOld:" << computeL2ErrorOld(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << " ----- L2NORM ----" << std::endl;
+ auto L2Norm = computeL2Norm(Basis_CE,phi_1);
+ std::cout << "L2Norm(phi_1): " << L2Norm << std::endl;
std::cout << " -----------------" << std::endl;
-
-
-
-
+
log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
-
-
+ log << "L2-Norm(Symphi_1): " << L2Norm_Symphi<< std::endl;
+ log << "L2-Norm(SymAnalytic): " << L2Norm_SymAnalytic << std::endl;
double EOC = 0.0;
@@ -2605,8 +1800,8 @@ int main(int argc, char *argv[])
// Storage_Error.push_back(5);
// Storage_Error.push_back("Second Entry");
- std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter-1)*6)+1 << std::endl; // Besser std::map ???
- std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl; // muss Storage_Error[idx] muss idx zur compile time feststehen?!
+ std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter-1)*6)+1 << std::endl; // Besser std::map ???
+ std::cout << " ((levelCounter-1)*6)+1: " << ((levelCounter)*6)+1 << std::endl; // für Storage_Error[idx] muss idx zur compile time feststehen?!
auto ErrorOld = std::get<double>(Storage_Error[((levelCounter-1)*6)+1]);
auto ErrorNew = std::get<double>(Storage_Error[(levelCounter*6)+1]);
@@ -2615,30 +1810,18 @@ int main(int argc, char *argv[])
// std::cout << "Storage_Error[0] :" << std::get<1>(Storage_Error[0]) << std::endl;
// std::cout << "output variant :" << std::get<std::string>(Storage_Error[1]) << std::endl;
// std::cout << "output variant :" << std::get<0>(Storage_Error[1]) << std::endl;
-
}
-
-
-
// Storage_Error.push_back(level);
Storage_Error.push_back(EOC);
Storage_Error.push_back(L2Norm_Symphi);
Storage_Error.push_back(L2Norm_SymAnalytic);
Storage_Error.push_back(L2Norm);
-// Storage_Error.push_back();
-// Storage_Error.push_back();
-// Storage_Error.push_back();
-// Storage_Error.push_back();
-
-
}
-
//////////////////////////////////////////////////////////////////////////////////////////////
// Write Data to Matlab / Optimization
//////////////////////////////////////////////////////////////////////////////////////////////
-
// writeMatrixToMatlab(Q, "matlab.txt");
writeMatrixToMatlab(Q, "../../Matlab-Programs/QMatrix.txt");
@@ -2648,19 +1831,13 @@ int main(int argc, char *argv[])
FieldMatrix<double,1,3> BeffMat;
BeffMat[0] = Beff;
-writeMatrixToMatlab(BeffMat, "../../Matlab-Programs/BMatrix.txt");
-
-
-
-
-
-
+ writeMatrixToMatlab(BeffMat, "../../Matlab-Programs/BMatrix.txt");
+
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
//////////////////////////////////////////////////////////////////////////////////////////////
-
VTKWriter<GridView> vtkWriter(gridView_CE);
vtkWriter.addVertexData(