diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 22e990f66eb4e1606c121382d67935e059840b8f..d52e73f334475c7f1ba39e77371ef1f827712937 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -49,7 +49,7 @@
#include <dune/functions/functionspacebases/hierarchicvectorwrapper.hh>
// #include <dune/fufem/discretizationerror.hh>
-// #include <boost/multiprecision/cpp_dec_float.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
using namespace Dune;
@@ -740,6 +740,10 @@ auto energy(const Basis& basis,
const MatrixFunction& matrixFieldFuncB)
{
+ // TEST HIGHER PRECISION
+// using float_50 = boost::multiprecision::cpp_dec_float_50;
+// float_50 higherPrecEnergy = 0.0;
+
double energy = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -772,8 +776,11 @@ auto energy(const Basis& basis,
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();
const auto& localFiniteElement = localView.tree().child(0).finiteElement();
@@ -785,7 +792,8 @@ auto energy(const Basis& basis,
// 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 double integrationElement = geometry.integrationElement(quadPos);
@@ -793,7 +801,7 @@ auto energy(const Basis& basis,
auto strain1 = matrixFieldA(quadPos);
// 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 ?
double trace = 0.0;
@@ -813,9 +821,13 @@ auto energy(const Basis& basis,
// elementEnergy += energyDensity * quad[pt].weight() * integrationElement;
elementEnergy += energyDensity * quadPoint.weight() * integrationElement;
+// elementEnergy_HP += energyDensity * quadPoint.weight() * integrationElement;
}
energy += elementEnergy;
+// higherPrecEnergy += elementEnergy_HP;
}
+ //TEST
+// std::cout << std::setprecision(std::numeric_limits<float_50>::digits10) << higherPrecEnergy << std::endl;
return energy;
}
@@ -1276,6 +1288,9 @@ double computeL2SymNormCoeff(const Basis& basis,
for (size_t k=0; k < nCompo; k++)
+
+
+
for (size_t i=0; i < nSf; i++)
{
for (size_t l=0; l < nCompo; l++)
@@ -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
template<class Basis, class Vector>
@@ -1339,6 +1453,157 @@ double computeL2NormCoeff(const Basis& basis,
constexpr int nCompo = 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 < 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 Domain = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
// using MatrixRT = FieldMatrix< double, nCompo, nCompo>;
@@ -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>
@@ -1616,7 +1930,7 @@ int main(int argc, char *argv[])
unsigned int prestraintype = parameterSet.get<unsigned int>("prestrainType", 2);
double rho1 = parameterSet.get<double>("rho1", 1.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;
auto prestrainImp = PrestrainImp(rho1, rho2, theta, width);
@@ -1676,7 +1990,9 @@ 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))
+// 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;
@@ -1686,7 +2002,7 @@ int main(int argc, char *argv[])
};
auto lambdaTerm = [lambda1,lambda2, theta] (const Domain& z) {
- if (abs(z[0]) >= (theta/2))
+ if (abs(z[0]) >= (theta/2.0))
return lambda1;
else
return lambda2;
@@ -1803,6 +2119,23 @@ int main(int argc, char *argv[])
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[])
// printvector(std::cout, row, "row" , "--");
//
-
+
@@ -2107,6 +2440,11 @@ int main(int argc, char *argv[])
double GGterm = 0.0;
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]*loadContainer[b]) + GGterm; // TEST
@@ -2216,20 +2554,25 @@ int main(int argc, char *argv[])
// 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
+ 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;
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;
+ auto L2Norm_SymAnalytic = computeL2SymError(Basis_CE,zeroVec ,gamma,symPhi_1_analytic);
+ std::cout << "L2-Norm(Symanalytic): " << L2Norm_SymAnalytic << std::endl;
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 << "L2Norm(phi_1) - GVfunc: " << computeL2NormFunc(Basis_CE,phi_1) << std::endl;
+ 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;
}
//////////////////////////////////////////////////////////////////////////////////////////////