diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 58d7d1deac04a0a835720ab87d61ed9ebe7e64dd..34e48b38d69cf419fd822d113cb579a337f444d4 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -1025,8 +1025,8 @@ double computeL2SymErrorNew(const Basis& basis,
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 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()))
@@ -1062,7 +1062,6 @@ double computeL2SymErrorNew(const Basis& basis,
MatrixRT tmp(0);
- MatrixRT tmpNew(0);
double sum = 0.0;
for (size_t k=0; k < nCompo; k++)
@@ -1082,44 +1081,21 @@ double computeL2SymErrorNew(const Basis& basis,
// std::cout << "coeffVector[globalIdx1]" << coeffVector[globalIdx1] << std::endl;
// 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]);
-// }
+ 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;
- tmpNew += defGradientU;
// 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;
-// }
- }
-
-// printmatrix(std::cout, tmpNew, "tmpNew", "--");
-
- for (int ii=0; ii<nCompo; ii++)
- for (int jj=0; jj<nCompo; jj++)
- {
- tmp[ii][jj] = 0.5 * (tmpNew[ii][jj] + tmpNew[jj][ii]);
+ tmp += strainU;
}
-
-
+
// printmatrix(std::cout, matrixField(quadPos), "matrixField(quadPos)", "--");
// printmatrix(std::cout, tmp, "tmp", "--"); // TEST Symphi_1 hat ganz andere Struktur als analytic?
-
// tmp = tmp - matrixField(quadPos);
// printmatrix(std::cout, tmp - matrixField(quadPos), "Difference", "--");
@@ -1128,15 +1104,12 @@ double computeL2SymErrorNew(const Basis& basis,
for (int jj=0; jj<nCompo; jj++)
{
sum += std::pow(tmp[ii][jj]-matrixField(quadPos)[ii][jj],2);
-// sum += std::pow(tmp[ii][jj],2);
}
// std::cout << "sum:" << sum << std::endl;
-
error += sum * quadPoint.weight() * integrationElement;
// std::cout << "error:" << error << std::endl;
}
}
-
return sqrt(error);
}
@@ -1153,6 +1126,7 @@ double computeL2SymError(const Basis& basis,
const MatrixFunction& matrixFieldFunc
)
{
+ // This Version computes the SCALAR PRODUCT of the difference ..
double error = 0.0;
constexpr int dim = 3;
constexpr int nCompo = 3;
@@ -1212,7 +1186,6 @@ double computeL2SymError(const Basis& basis,
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?
@@ -1252,17 +1225,17 @@ 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
+ 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<nCompo; ii++)
for (int jj=0; jj<nCompo; jj++)
{
- strainU[ii][jj] = 0.5 * (defGradientU_2[ii][jj] + defGradientU_2[jj][ii]);
+ strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
}
double tmp = 0.0;
@@ -2040,6 +2013,105 @@ 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 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>;
+
+ 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 < nCompo; 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<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ {
+ strainU[ii][jj] = 0.5 * (defGradientU[ii][jj] + defGradientU[jj][ii]);
+ }
+
+ double tmp = 0.0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ tmp += (-2.0)*coeffVector[globalIdx1]*strainU[ii][jj] * matrixField(quadPos)[ii][jj];
+
+ error += tmp * quadPoint.weight() * integrationElement;
+ }
+
+ }
+ }
+
+
+ return sqrt(error+Term1+Term2);
+}
+
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -2726,6 +2798,7 @@ int main(int argc, char *argv[])
std::cout << " ----- TEST ----" << std::endl;
std::cout << "computeL2SymErrorNew: " << computeL2SymErrorNew(Basis_CE,phi_1,gamma,symPhi_1_analytic) << std::endl;
+ std::cout << "computeL2SymErrorSPTest: " << computeL2SymErrorSPTest(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);