From 1e2e0adcae52a359e4dea80d6ad4b33c44c6eb28 Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Thu, 5 Aug 2021 12:04:30 +0200
Subject: [PATCH] Add L2-Norm & run checks
---
src/dune-microstructure.cc | 190 +++++++++++++++++++++++++++++++++++--
1 file changed, 181 insertions(+), 9 deletions(-)
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 9fee3c16..22e990f6 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -49,6 +49,7 @@
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
// #include <dune/fufem/discretizationerror.hh>
+// #include <boost/multiprecision/cpp_dec_float.hpp>
using namespace Dune;
@@ -793,7 +794,7 @@ auto energy(const Basis& basis,
// St.Venant-Kirchhoff stress
MatrixRT stressU(0);
- stressU.axpy(2*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
+ stressU.axpy(2.0*muLocal(quadPos), strain1); //this += 2mu *strainU // eigentlich this += 2mu *strain1 ?
double trace = 0.0;
for (int ii=0; ii<nCompo; ii++)
@@ -1182,6 +1183,7 @@ double computeL2SymError(const Basis& basis,
}
}
+ MatrixRT defGradientU_2(0); //TEST (doesnt change anything..)
for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++)
{
@@ -1189,16 +1191,16 @@ double computeL2SymError(const Basis& basis,
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
+ defGradientU_2[k][0] = gradients[i][0]; // Y
+ defGradientU_2[k][1] = gradients[i][1]; //X2
+ defGradientU_2[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_2[ii][jj] + defGradientU_2[jj][ii]);
}
double tmp = 0.0;
@@ -1223,9 +1225,175 @@ double computeL2SymError(const Basis& basis,
+// TEST
+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 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();
+// 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 < nCompo; k++)
+ for (size_t i=0; i < nSf; i++)
+ {
+ 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); // 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
+ 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]);
+ }
+
+ // 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];
+
+
+ error += tmp * quadPoint.weight() * integrationElement;
+ }
+ }
+
+
+ }
+ }
+
+ return sqrt(error);
+}
+
+
+
+// TODO
+
+
+// TEST
+template<class Basis, class Vector>
+double computeL2NormCoeff(const Basis& basis,
+ Vector& coeffVector
+ )
+{
+ double error = 0.0;
+ 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();
+// 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 < nCompo; 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);
+}
@@ -1548,7 +1716,7 @@ int main(int argc, char *argv[])
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
/////////////////////////////////////////////////////////
@@ -1733,6 +1901,7 @@ int main(int argc, char *argv[])
solver.apply(x_2, load_alpha2, statistics);
std::cout << "solve linear system for x_3.\n";
solver.apply(x_3, load_alpha3, statistics);
+ log << "Solver-type used: " <<"\n CG-Solver" << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
@@ -1763,6 +1932,7 @@ int main(int argc, char *argv[])
solver.apply(x_1, load_alpha1, statistics);
solver.apply(x_2, load_alpha2, statistics);
solver.apply(x_3, load_alpha3, statistics);
+ log << "Solver-type used: " <<"\n GMRES-Solver" << std::endl;
}
////////////////////////////////////////////////////////////////////////////////////
else if ( Solvertype == 3)// QR - SOLVER
@@ -1778,7 +1948,7 @@ int main(int argc, char *argv[])
sPQR.apply(x_1, load_alpha1, statisticsQR);
sPQR.apply(x_2, load_alpha2, statisticsQR);
sPQR.apply(x_3, load_alpha3, statisticsQR);
-
+ log << "Solver-type used: " <<"\n QR-Solver" << std::endl;
}
// printvector(std::cout, load_alpha1BS, "load_alpha1 before SOLVER", "--" );
// printvector(std::cout, load_alpha1, "load_alpha1 AFTER SOLVER", "--" );
@@ -2009,7 +2179,7 @@ int main(int argc, char *argv[])
double b1 = (-(theta/4.0)*mu1*mu2)/(theta*mu1+(1.0-theta)*mu2);
double b2 = -(theta/4.0)*mu2;
double b3 = 0.0;
- double q1 = (mu1/6.0)*mu2/(theta*mu1+ (1.0- theta)*mu2);
+ double q1 = ((mu1*mu2)/6.0)/(theta*mu1+ (1.0- theta)*mu2);
double q2 = ((1.0-theta)*mu1+theta*mu2)/6.0;
std::cout << " --- print analytic solutions --- " << std::endl;
@@ -2029,7 +2199,7 @@ int main(int argc, char *argv[])
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}};
+ 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) {
@@ -2058,6 +2228,8 @@ int main(int argc, char *argv[])
log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
//TEST
+ std::cout << "L2Norm(phi_1):" << computeL2NormCoeff(Basis_CE,phi_1) << std::endl;
+ std::cout << "L2Norm(SymPhi_1): " << computeL2SymNormCoeff(Basis_CE,phi_1,gamma) << std::endl; // does not make a difference
std::cout << "L2ErrorOld:" << computeL2ErrorOld(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
}
//////////////////////////////////////////////////////////////////////////////////////////////
--
GitLab