diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index d52e73f334475c7f1ba39e77371ef1f827712937..31a59ba66eb6e414ed4ea52cb6dde193d06d5d04 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -1097,6 +1097,118 @@ double evalSymGrad(const Basis& basis,
+template<class Basis, class Vector, class MatrixFunction>
+double computeL2SymErrorNew(const Basis& basis,
+ Vector& coeffVector,
+ const double gamma,
+ const MatrixFunction& matrixFieldFunc
+ )
+{
+ double error = 0.0;
+ constexpr int dim = 3;
+ constexpr int nCompo = 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, nCompo, nCompo>;
+
+ 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]);
+
+
+ MatrixRT tmp(0);
+
+ double sum = 0.0;
+
+ for (size_t i=0; i < nSf; i++)
+ for (size_t k=0; k < nCompo; 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<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, tmp, "tmp", "--");
+
+
+ // 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[ii][jj] += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
+//
+//
+// error += tmp * quadPoint.weight() * integrationElement;
+// }
+ }
+
+ tmp = tmp - matrixField(quadPos);
+
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ sum += std::pow(tmp[ii][jj],2);
+
+ error += sum * quadPoint.weight() * integrationElement;
+
+ }
+ }
+
+ return sqrt(error);
+}
+
+
+
+
+
template<class Basis, class Vector, class MatrixFunction>
@@ -1134,7 +1246,6 @@ double computeL2SymError(const Basis& basis,
for (const auto& quadPoint : quad)
{
-
const auto& quadPos = quadPoint.position();
const auto jacobianInverseTransposed = geometry.jacobianInverseTransposed(quadPoint.position());
const auto integrationElement = geometry.integrationElement(quadPoint.position());
@@ -1148,10 +1259,10 @@ double computeL2SymError(const Basis& basis,
// std::vector< VectorRT> gradients(referenceGradients.size());
for (size_t i=0; i<gradients.size(); i++)
- jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
+ jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]); // TRANSFORMED
- MatrixRT defGradientU(0), defGradientV(0);
+// MatrixRT defGradientU(0), defGradientV(0);
for (size_t k=0; k < nCompo; k++)
@@ -1165,12 +1276,15 @@ double computeL2SymError(const Basis& basis,
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
@@ -1195,7 +1309,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++)
{
@@ -1203,6 +1317,7 @@ double computeL2SymError(const Basis& basis,
size_t globalIdx1 = localView.index(localIdx1);
// (scaled) Deformation gradient of the ansatz basis function
+ MatrixRT defGradientU_2(0); //TEST (doesnt change anything..)
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
@@ -1237,6 +1352,8 @@ double computeL2SymError(const Basis& basis,
+
+
// TEST
template<class Basis, class Vector>
double computeL2SymNormCoeff(const Basis& basis,
@@ -1284,7 +1401,116 @@ double computeL2SymNormCoeff(const Basis& basis,
jacobianInverseTransposed.mv(referenceGradients[i][0], gradients[i]);
- MatrixRT defGradientU(0), defGradientV(0);
+ MatrixRT tmp(0);
+
+ double sum = 0.0;
+
+ for (size_t i=0; i < nSf; i++)
+ for (size_t k=0; k < nCompo; 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<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, tmp, "tmp", "--");
+
+
+ // 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[ii][jj] += coeffVector[globalIdx1]*coeffVector[globalIdx2]* strainU[ii][jj] * strainV[ii][jj];
+//
+//
+// error += tmp * quadPoint.weight() * integrationElement;
+// }
+ }
+
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; 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 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++)
@@ -1303,10 +1529,12 @@ double computeL2SymNormCoeff(const Basis& basis,
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
@@ -1345,7 +1573,7 @@ double computeL2SymNormCoeff(const Basis& basis,
template<class Basis, class Vector>
-double computeL2SymNormCoeff(const Basis& basis,
+double computeL2SymNormCoeffV3(const Basis& basis,
Vector& coeffVector,
const double gamma
)
@@ -1442,6 +1670,8 @@ double computeL2SymNormCoeff(const Basis& basis,
}
*/
+
+ ////////////////////////////////////////////////////////////// L2-NORM /////////////////////////////////////////////////////////////////
// TEST
template<class Basis, class Vector>
double computeL2NormCoeff(const Basis& basis,
@@ -2547,6 +2777,12 @@ int main(int argc, char *argv[])
if(write_L2Error)
{
+
+ std::cout << " ----- TEST ----" << std::endl;
+
+ std::cout << "computeL2SymErrorNew: " << computeL2SymError(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);
std::cout << "L2-SymError: " << L2SymError << std::endl;
@@ -2556,6 +2792,7 @@ int main(int argc, char *argv[])
auto L2Norm_Symphi = computeL2SymError(Basis_CE,phi_1,gamma,zeroFunction); // Achtung Norm SymPhi_1
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;