From 862e8d488afcc4e112cd57505d6ba9b2e7432f21 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20H=C3=B6rl?= <maximilian.hoerl@mathematik.uni-stuttgart.de> Date: Thu, 13 Feb 2020 19:20:21 +0100 Subject: [PATCH] use quadrature rule for quadratic integrals in dg.hh --- dune/mmdg/dg.hh | 57 +++++++++++++++++++++++-------------------------- 1 file changed, 27 insertions(+), 30 deletions(-) diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh index 76def0c..fb2e8c7 100644 --- a/dune/mmdg/dg.hh +++ b/dune/mmdg/dg.hh @@ -78,6 +78,7 @@ private: QRulesEdge::rule(edge, problem_.quadratureOrder_K()); const QRuleEdge& rule_Kplus1_Edge = QRulesEdge::rule(edge, problem_.quadratureOrder_K() + 1); + const QRuleEdge& rule_2ndOrder = QRulesEdge::rule(edge, 2); //we use the basis // phi_elem,0 (x) = indicator(elem); @@ -147,11 +148,6 @@ private: const double intersctVol = intersctGeo.volume(); const auto& intersctCenter = intersctGeo.center(); - //TODO: quadrature rule cannot be used for dim = 1! - // const Dune::QuadratureRule<double,dim-1>& secondOrderRule = - // Dune::QuadratureRules<double,dim-1>::rule( - // intersct.type(), 2); - //create storage for multiply used integral values //linearIntegrals[i] = int_intersct x[i] ds Dune::FieldVector<Scalar, dim> linearIntegrals(0.0); @@ -168,37 +164,38 @@ private: //KIntegralsX[i][j] = int_intersct x[j]*(K[:,i]*normal) ds Dune::FieldMatrix<Scalar, dim, dim> KIntegralsX(0.0); - //fill vector linearIntegrals and matrix quadraticIntregrals //NOTE - for (int i = 0; i < dim; i++) + //fill vector linearIntegrals and matrix quadraticIntregrals + for (const auto& ip : rule_2ndOrder) { - //use midpoint rule for exact evaluation of - // int_intersct x[i] ds - linearIntegrals[i] = intersctVol * intersctCenter[i]; + //quadrature point in local and global coordinates + const auto& qpLocal = ip.position(); + const auto& qpGlobal = intersctGeo.global(qpLocal); - for (int j = 0; j <= i; j++) + const Scalar weight(ip.weight()); + const Scalar integrationElem(intersctGeo.integrationElement(qpLocal)); + + for (int i = 0; i < dim; i++) { - const auto& leftCorner = intersctGeo.corner(0); - const auto& rightCorner = intersctGeo.corner(1); - - quadraticIntregrals[i][j] = intersctVol / 3 * - ( leftCorner[i] * leftCorner[j] + rightCorner[i] * rightCorner[j] - + 0.5 * (leftCorner[i] * rightCorner[j] + - leftCorner[j] * rightCorner[i]) ); - /* - //use second order quadrature rule for exact evaluation of - // int_intersct x_i*x_j ds - //NOTE: use Simpson's rule instead manually? - for (const auto& ip : secondOrderRule) + //use midpoint rule for exact evaluation of + // int_intersct x[i] ds + linearIntegrals[i] = intersctVol * intersctCenter[i]; + + for (int j = 0; j <= i; j++) { - const auto& qp = ip.position(); - quadraticIntregrals[i][j] += //NOTE: volume necessary? - ip.weight() * qp[i] * qp[j] * intersctVol; + quadraticIntregrals[i][j] += + weight * integrationElem * qpGlobal[i] * qpGlobal[j]; + } + } + } + + //quadraticIntregrals is symmetric + for (int i = 0; i < dim; i++) + { + for (int j = 0; j < i; j++) + { + quadraticIntregrals[j][i] = quadraticIntregrals[i][j]; } - */ - //NOTE: probably unnecessary - quadraticIntregrals[j][i] = quadraticIntregrals[i][j]; } - } //fill vector KIntegrals using a quadrature rule for (const auto& ip : rule_K_Edge) -- GitLab