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

added L2-Error

parent 5e38dc69
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/dune-microstructure.cc
+ 137
− 8
View file @ 029bcc6e
......@@ -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
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment