Skip to content
Snippets Groups Projects
Commit 0af9b0de authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

3 different Versions of L2-Norm

parent 690aadba
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/dune-microstructure.cc
+ 358
− 15
View file @ 0af9b0de
...@@ -49,7 +49,7 @@ ...@@ -49,7 +49,7 @@
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh> #include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
// #include <dune/fufem/discretizationerror.hh> // #include <dune/fufem/discretizationerror.hh>
// #include <boost/multiprecision/cpp_dec_float.hpp> #include <boost/multiprecision/cpp_dec_float.hpp>
using namespace Dune; using namespace Dune;
...@@ -740,6 +740,10 @@ auto energy(const Basis& basis, ...@@ -740,6 +740,10 @@ auto energy(const Basis& basis,
const MatrixFunction& matrixFieldFuncB) const MatrixFunction& matrixFieldFuncB)
{ {
// TEST HIGHER PRECISION
// using float_50 = boost::multiprecision::cpp_dec_float_50;
// float_50 higherPrecEnergy = 0.0;
double energy = 0.0; double energy = 0.0;
constexpr int dim = 3; constexpr int dim = 3;
constexpr int nCompo = 3; constexpr int nCompo = 3;
...@@ -772,8 +776,11 @@ auto energy(const Basis& basis, ...@@ -772,8 +776,11 @@ auto energy(const Basis& basis,
lambdaLocal.bind(e); lambdaLocal.bind(e);
FieldVector<double,1> elementEnergy(0); // FieldVector<double,1> elementEnergy(0); //!!!
// printvector(std::cout, elementEnergy, "elementEnergy" , "--");
double elementEnergy = 0.0;
// double elementEnergy_HP = 0.0;
auto geometry = e.geometry(); auto geometry = e.geometry();
const auto& localFiniteElement = localView.tree().child(0).finiteElement(); const auto& localFiniteElement = localView.tree().child(0).finiteElement();
...@@ -785,7 +792,8 @@ auto energy(const Basis& basis, ...@@ -785,7 +792,8 @@ auto energy(const Basis& basis,
// for (size_t pt=0; pt < quad.size(); pt++) { // for (size_t pt=0; pt < quad.size(); pt++) {
const FieldVector<double,dim>& quadPos = quadPoint.position(); // const FieldVector<double,dim>& quadPos = quadPoint.position();
const auto& quadPos = quadPoint.position();
// const Domain& quadPos = quad[pt].position(); // const Domain& quadPos = quad[pt].position();
const double integrationElement = geometry.integrationElement(quadPos); const double integrationElement = geometry.integrationElement(quadPos);
...@@ -793,7 +801,7 @@ auto energy(const Basis& basis, ...@@ -793,7 +801,7 @@ auto energy(const Basis& basis,
auto strain1 = matrixFieldA(quadPos); auto strain1 = matrixFieldA(quadPos);
// St.Venant-Kirchhoff stress // St.Venant-Kirchhoff stress
MatrixRT stressU(0); MatrixRT stressU(0.0);
stressU.axpy(2.0*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; double trace = 0.0;
...@@ -813,9 +821,13 @@ auto energy(const Basis& basis, ...@@ -813,9 +821,13 @@ auto energy(const Basis& basis,
// elementEnergy += energyDensity * quad[pt].weight() * integrationElement; // elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
elementEnergy += energyDensity * quadPoint.weight() * integrationElement; elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
// elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
} }
energy += elementEnergy; energy += elementEnergy;
// higherPrecEnergy += elementEnergy_HP;
} }
//TEST
// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
return energy; return energy;
} }
...@@ -1276,6 +1288,9 @@ double computeL2SymNormCoeff(const Basis& basis, ...@@ -1276,6 +1288,9 @@ double computeL2SymNormCoeff(const Basis& basis,
for (size_t k=0; k < nCompo; k++) for (size_t k=0; k < nCompo; k++)
for (size_t i=0; i < nSf; i++) for (size_t i=0; i < nSf; i++)
{ {
for (size_t l=0; l < nCompo; l++) for (size_t l=0; l < nCompo; l++)
...@@ -1324,9 +1339,108 @@ double computeL2SymNormCoeff(const Basis& basis, ...@@ -1324,9 +1339,108 @@ double computeL2SymNormCoeff(const Basis& basis,
} }
/*
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;
}
}
// TODO
}
}
return sqrt(error);
}
*/
// TEST // TEST
template<class Basis, class Vector> template<class Basis, class Vector>
...@@ -1339,6 +1453,157 @@ double computeL2NormCoeff(const Basis& basis, ...@@ -1339,6 +1453,157 @@ double computeL2NormCoeff(const Basis& basis,
constexpr int nCompo = 3; constexpr int nCompo = 3;
auto localView = basis.localView(); 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 < nCompo; 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 < 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(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 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++) // ANALOG STOKES Beitrag nur wenn k=l!!!
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(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 nCompo = 3;
auto localView = basis.localView();
// using GridView = typename Basis::GridView; // using GridView = typename Basis::GridView;
// using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate; // using Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
// using MatrixRT = FieldMatrix< double, nCompo, nCompo>; // using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
...@@ -1396,6 +1661,55 @@ double computeL2NormCoeff(const Basis& basis, ...@@ -1396,6 +1661,55 @@ double computeL2NormCoeff(const Basis& basis,
} }
// TEST
template<class Basis, class Vector>
double computeL2NormFunc(const Basis& basis,
Vector& coeffVector
)
{
// IMPLEMENTATION with makeDiscreteGlobalBasisFunction
double error = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
auto localView = basis.localView();
using SolutionRange = FieldVector<double,dim>;
auto GVFunc = Functions::makeDiscreteGlobalBasisFunction<SolutionRange>(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);
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);
for(const auto& quadPoint : quad)
{
const double integrationElement = element.geometry().integrationElement(quadPoint.position());
const auto& quadPos = quadPoint.position();
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> template<class Basis, class Vector, class MatrixFunction>
...@@ -1616,7 +1930,7 @@ int main(int argc, char *argv[]) ...@@ -1616,7 +1930,7 @@ int main(int argc, char *argv[])
unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", 2); unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", 2);
double rho1 = parameterSet.get<double>("rho1", 1.0); double rho1 = parameterSet.get<double>("rho1", 1.0);
double alpha = parameterSet.get<double>("alpha", 2.0); double alpha = parameterSet.get<double>("alpha", 2.0);
double theta = parameterSet.get<double>("theta",0.3); double theta = parameterSet.get<double>("theta",0.3); //TEST theta = 1.0/4.0
double rho2 = alpha*1.0; double rho2 = alpha*1.0;
auto prestrainImp = PrestrainImp(rho1, rho2, theta, width); auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
...@@ -1676,7 +1990,9 @@ int main(int argc, char *argv[]) ...@@ -1676,7 +1990,9 @@ int main(int argc, char *argv[])
// Create Lambda-Functions for material Parameters depending on microstructure // Create Lambda-Functions for material Parameters depending on microstructure
/////////////////////////////////// ///////////////////////////////////
auto muTerm = [mu1, mu2, theta] (const Domain& z) { auto muTerm = [mu1, mu2, theta] (const Domain& z) {
if (abs(z[0]) >= (theta/2)) // 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; return mu1;
// if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das? // if (abs(z[0]) < (theta/2) && z[2] < 0) // oder das?
// return mu2; // return mu2;
...@@ -1686,7 +2002,7 @@ int main(int argc, char *argv[]) ...@@ -1686,7 +2002,7 @@ int main(int argc, char *argv[])
}; };
auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) { auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
if (abs(z[0]) >= (theta/2)) if (abs(z[0]) >= (theta/2.0))
return lambda1; return lambda1;
else else
return lambda2; return lambda2;
...@@ -1803,6 +2119,23 @@ int main(int argc, char *argv[]) ...@@ -1803,6 +2119,23 @@ int main(int argc, char *argv[])
return -1.0*x3G_3(x); return -1.0*x3G_3(x);
}; };
// TEST - energy method ///
// different indicatorFunction in muTerm has impact here !!
// double GGterm = 0.0;
// GGterm = energy(Basis_CE, muLocal, lambdaLocal, x3G_1, x3G_1 ); // <L i(x3G_alpha) , i(x3G_beta) >
// std::cout << "GGTerm:" << GGterm << std::endl;
// std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
//
//////////////////////////////////////////////////
/////////////////////////////////////////////// ///////////////////////////////////////////////
...@@ -1852,7 +2185,7 @@ int main(int argc, char *argv[]) ...@@ -1852,7 +2185,7 @@ int main(int argc, char *argv[])
// printvector(std::cout, row, "row" , "--"); // printvector(std::cout, row, "row" , "--");
// //
...@@ -2107,6 +2440,11 @@ int main(int argc, char *argv[]) ...@@ -2107,6 +2440,11 @@ int main(int argc, char *argv[])
double GGterm = 0.0; 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) >
// TEST
//TEST
std::cout << " analyticGGTERM:" << (mu1*(1-theta)+mu2*theta)/6.0 << std::endl;
std::setprecision(std::numeric_limits<float>::digits10);
Q[a][b] = (coeffContainer[a]*tmp1) + GGterm; // seems symmetric... check positiv definitness? 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
...@@ -2216,20 +2554,25 @@ int main(int argc, char *argv[]) ...@@ -2216,20 +2554,25 @@ int main(int argc, char *argv[])
// std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl; // std::cout << "L2-ErrorTEST: " << L2errorTest << std::endl;
auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1 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 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
VectorCT zeroVec; VectorCT zeroVec;
zeroVec.resize(Basis_CE.size()); zeroVec.resize(Basis_CE.size());
zeroVec = 0; zeroVec = 0;
auto L2Norm_analytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic); auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
std::cout << "L2-Norm(analytic): " << L2Norm_analytic << std::endl; std::cout << "L2-Norm(Symanalytic): " << L2Norm_SymAnalytic << std::endl;
log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl; log << "L2-Error (symmetric Gradient phi_1):" << L2SymError << std::endl;
//TEST //TEST
std::cout << "L2Norm(phi_1):" << computeL2NormCoeff(Basis_CE,phi_1) << std::endl; std::cout << "L2Norm(phi_1) - GVfunc: " << computeL2NormFunc(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 << "L2Norm(phi_1) - coeff: " << computeL2NormCoeff(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 << "L2ErrorOld:" << computeL2ErrorOld(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
} }
////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment