From 2bed543146cfe6e8de50700706065cd8271117ac Mon Sep 17 00:00:00 2001
From: Klaus <klaus.boehnlein@tu-dresden.de>
Date: Thu, 15 Jul 2021 20:46:54 +0200
Subject: [PATCH] added L2-Error
---
src/dune-microstructure.cc | 145 +++++++++++++++++++++++++++++++++++--
1 file changed, 137 insertions(+), 8 deletions(-)
diff --git a/src/dune-microstructure.cc b/src/dune-microstructure.cc
index 877b8136..5669c14a 100644
--- a/src/dune-microstructure.cc
+++ b/src/dune-microstructure.cc
@@ -106,12 +106,12 @@ void computeElementStiffnessMatrix(const LocalView& localView,
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
//print:
// printmatrix(std::cout, basisContainer[0] , "G_1", "--");
// printmatrix(std::cout, basisContainer[1] , "G_2", "--");
- // printmatrix(std::cout, basisContainer[2] , "G_3", "--");
+// printmatrix(std::cout, basisContainer[2] , "G_3", "--");
////////////////////////////////////////////////////
int orderQR = 2*(dim*localFiniteElement.localBasis().order()-1);
@@ -501,7 +501,7 @@ void computeElementLoadVector( const LocalView& localView,
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
@@ -1220,11 +1220,105 @@ auto equivalent = [](const FieldVector<double,3>& x, const FieldVector<double,3>
+template<class Basis, class Vector, class MatrixFunction>
+double computeL2Error(const Basis& basis,
+ Vector& coeffVector,
+ const double gamma,
+ const MatrixFunction& matrixFieldFunc
+ )
+{
+
+ auto 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();
+
+
+ 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);
+
+
+ for (size_t k=0; k < nCompo; k++)
+ for (size_t i=0; i < nSf; i++)
+ {
+ // (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
+ }
+
+
+ // 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]); // symmetric gradient
+ }
+
+
+ // Local energy density: stress times strain
+ double tmp = 0;
+ for (int ii=0; ii<nCompo; ii++)
+ for (int jj=0; jj<nCompo; jj++)
+ tmp+= std::pow(strainU[ii][jj] - matrixField(quadPos)[ii][jj] ,2);
+
+
+
+ error += tmp * quadPoint.weight() * integrationElement;
+ }
+ }
+
+ return sqrt(error);
+
+}
@@ -1350,8 +1444,10 @@ int main(int argc, char *argv[])
+// double beta = 1;
double beta = 2;
-
+
+
double mu1 = 10;
// double mu1 = 0.5*17e6;
@@ -1482,10 +1578,12 @@ int main(int argc, char *argv[])
};
Func2Tensor x3G_3 = [] (const Domain& x) {
- return MatrixRT{{0.0, 1/sqrt(2)*x[2], 0.0}, {1/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ return MatrixRT{{0.0, 1.0/sqrt(2)*x[2], 0.0}, {1.0/sqrt(2)*x[2], 0.0, 0.0}, {0.0, 0.0, 0.0}};
};
-
+// FieldVector<double,3> test= {1.0/4.0 , 0.0 , 0.25};
+// auto Tast = x3G_3(test);
+// printmatrix(std::cout, Tast, "x3G_3", "--");
/////////////////////////////////////////////////////////////////////// TODO
// TODO : PrestrainImp.hh
@@ -1518,7 +1616,7 @@ int main(int argc, char *argv[])
//////////////////////////////////////////////
MatrixRT G_1 {{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0, 0.0}};
MatrixRT G_2 {{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0, 0.0, 0.0}};
- MatrixRT G_3 {{0.0, 1/sqrt(2), 0.0}, {1/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
+ MatrixRT G_3 {{0.0, 1.0/sqrt(2), 0.0}, {1.0/sqrt(2), 0.0, 0.0}, {0.0, 0.0, 0.0}};
std::array<MatrixRT,3 > basisContainer = {G_1, G_2, G_3};
@@ -1855,6 +1953,10 @@ int main(int argc, char *argv[])
log << "b2_hat: " << B_hat[2] << std::endl;
log << "b3_hat: " << B_hat[3] << std::endl;
+
+
+
+
//////////////////////////////////////////////////////////////
// Define Analytic Solutions
//////////////////////////////////////////////////////////////
@@ -1876,7 +1978,34 @@ int main(int argc, char *argv[])
std::cout << "q2 : " << q2 << std::endl;
std::cout << "q3 should be between q1 and q2" << std::endl;
-
+
+
+
+ // TODO Define sym grad phi_1
+ // - how to compute <mu>_h ?
+
+// Func2Tensor symPhi_1_analytic = [] (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}};
+// };
+
+ 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}};
+ };
+
+
+
+ FieldVector<double,3> testvector = {1.0/4.0 , 0.0 , 0.25};
+ std::cout << "t[2]: " << testvector[2] << std::endl;
+ std::cout << "muTerm(t):" << muTerm(testvector) << std::endl;
+ auto Teest = symPhi_1_analytic(testvector);
+ printmatrix(std::cout, Teest, "symPhi_1_analytic(t)", "--");
+
+
+
+
+ auto L2error = computeL2Error(Basis_CE,phi_1,gamma,symPhi_1_analytic);
+ std::cout << "L2-Error: " << L2error << std::endl;
+
//////////////////////////////////////////////////////////////////////////////////////////////
// Write result to VTK file
--
GitLab