From ba7a96565b1ebb15a4216ac6addedb3c41a7492f Mon Sep 17 00:00:00 2001 From: Nicola Giuliani Date: Fri, 17 Jan 2020 15:48:32 +0100 Subject: [PATCH] added numerical quadrature for non planar cells --- source/grid/tria_accessor.cc | 146 ++++++++++++++---- .../measure_of_distorted_codim_one_face.cc | 16 +- 2 files changed, 125 insertions(+), 37 deletions(-) diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index a7edc8283d..09015c006a 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1285,41 +1285,125 @@ namespace } + // a 2d face in 3d space + // a 2d face in 3d space template double - measure(const dealii::TriaAccessor<2, dim, 3> &accessor) + measure(const TriaAccessor<2, dim, 3> &accessor) { - // In general the area can be computed as - // the integral of the cross product of the two tangential vectors - - // If we assume a bilinear patch parametrized in u and v we get that - // t_u = (v_1 - v_0) + v (v_3 - v_2 - v_1 + v_0) - // t_v = (v_2 - v_0) + u (v_3 - v_2 - v_1 + v_0) - // So t_u x t_v = (v_1 - v_0) x (v_2 - v_0) + u (v_1 - v_0) x (v_3 - v_2 - - // v_1 + v_0) + v (v_3 - v_2 - v_1 + v_0) x (v_2 - v_0) t_u x t_v = w_1 + u - // w_2 + v w_3 we can integrate the square norm (t_u x t_v) * (t_u x t_v) = - // w_1*w_1 + u^2 w_2*w_2 + v^2 w_3*w_3 + 2u w_1*w_2 + 2v w_1*w_3 + 2uv - // w_2*w_3 in u and v getting (between zero and one) w_1*w_1 + 1/3 w_2*w_2 + - // 1/3 w_3*w_3 + w_1*w_2 + w_1*w_3 + 1/2 w_2*w_3 - - const Tensor<1, 3> w_1 = - cross_product_3d(accessor.vertex(1) - accessor.vertex(0), - accessor.vertex(2) - accessor.vertex(0)); - const Tensor<1, 3> w_2 = - cross_product_3d(accessor.vertex(1) - accessor.vertex(0), - accessor.vertex(3) - accessor.vertex(2) - - accessor.vertex(1) + accessor.vertex(0)); - const Tensor<1, 3> w_3 = - cross_product_3d(accessor.vertex(3) - accessor.vertex(2) - - accessor.vertex(1) + accessor.vertex(0), - accessor.vertex(2) - accessor.vertex(0)); - - - return std::sqrt(scalar_product(w_1, w_1) + scalar_product(w_1, w_2) + - scalar_product(w_1, w_3) + 0.5 * scalar_product(w_2, w_3) + - 1. / 3 * scalar_product(w_2, w_2) + - 1. / 3 * scalar_product(w_3, w_3)); + // If the face is planar, the diagonal from vertex 0 to vertex 3, + // v_03, should be in the plane P_012 of vertices 0, 1 and 2. Get + // the normal vector of P_012 and test if v_03 is orthogonal to + // that. If so, the face is planar and computing its area is simple. + const Tensor<1, 3> v01 = accessor.vertex(1) - accessor.vertex(0); + const Tensor<1, 3> v02 = accessor.vertex(2) - accessor.vertex(0); + + const Tensor<1, 3> normal = cross_product_3d(v01, v02); + + const Tensor<1, 3> v03 = accessor.vertex(3) - accessor.vertex(0); + + // check whether v03 does not lie in the plane of v01 and v02 + // (i.e., whether the face is not planar). we do so by checking + // whether the triple product (v01 x v02) * v03 forms a positive + // volume relative to |v01|*|v02|*|v03|. the test checks the + // squares of these to avoid taking norms/square roots: + if (std::abs((v03 * normal) * (v03 * normal) / + ((v03 * v03) * (v01 * v01) * (v02 * v02))) >= 1e-24) + { + // If the vectors are non planar we integrate the norm of the normal + // vector using a numerical Gauss scheme of order 4. + const Tensor<1, 3> w_1 = + cross_product_3d(accessor.vertex(1) - accessor.vertex(0), + accessor.vertex(2) - accessor.vertex(0)); + const Tensor<1, 3> w_2 = + cross_product_3d(accessor.vertex(1) - accessor.vertex(0), + accessor.vertex(3) - accessor.vertex(2) - + accessor.vertex(1) + accessor.vertex(0)); + const Tensor<1, 3> w_3 = + cross_product_3d(accessor.vertex(3) - accessor.vertex(2) - + accessor.vertex(1) + accessor.vertex(0), + accessor.vertex(2) - accessor.vertex(0)); + + double a = scalar_product(w_1, w_1); + double b = scalar_product(w_2, w_2); + double c = scalar_product(w_3, w_3); + double d = scalar_product(w_1, w_2); + double e = scalar_product(w_1, w_3); + double f = scalar_product(w_2, w_3); + + return 0.03025074832140047 * + std::sqrt(a + 0.0048207809894260144 * b + + 0.0048207809894260144 * c + 0.13886368840594743 * d + + 0.13886368840594743 * e + + 0.0096415619788520288 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.0048207809894260144 * b + + 0.10890625570683385 * c + 0.13886368840594743 * d + + 0.66001895641514374 * e + 0.045826333352825557 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.0048207809894260144 * b + + 0.44888729929169013 * c + 0.13886368840594743 * d + + 1.3399810435848563 * e + 0.09303735505312187 * f) + + 0.03025074832140047 * + std::sqrt(a + 0.0048207809894260144 * b + + 0.86595709258347853 * c + 0.13886368840594743 * d + + 1.8611363115940525 * e + 0.12922212642709538 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.10890625570683385 * b + + 0.0048207809894260144 * c + 0.66001895641514374 * d + + 0.13886368840594743 * e + 0.045826333352825557 * f) + + 0.10632332575267359 * + std::sqrt(a + 0.10890625570683385 * b + + 0.10890625570683385 * c + 0.66001895641514374 * d + + 0.66001895641514374 * e + 0.2178125114136677 * f) + + 0.10632332575267359 * + std::sqrt(a + 0.10890625570683385 * b + + 0.44888729929169013 * c + 0.66001895641514374 * d + + 1.3399810435848563 * e + 0.44220644500147605 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.10890625570683385 * b + + 0.86595709258347853 * c + 0.66001895641514374 * d + + 1.8611363115940525 * e + 0.61419262306231814 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.44888729929169013 * b + + 0.0048207809894260144 * c + 1.3399810435848563 * d + + 0.13886368840594743 * e + 0.09303735505312187 * f) + + 0.10632332575267359 * + std::sqrt(a + 0.44888729929169013 * b + + 0.10890625570683385 * c + 1.3399810435848563 * d + + 0.66001895641514374 * e + 0.44220644500147605 * f) + + 0.10632332575267359 * + std::sqrt(a + 0.44888729929169013 * b + + 0.44888729929169013 * c + 1.3399810435848563 * d + + 1.3399810435848563 * e + 0.89777459858338027 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.44888729929169013 * b + + 0.86595709258347853 * c + 1.3399810435848563 * d + + 1.8611363115940525 * e + 1.2469436885317342 * f) + + 0.03025074832140047 * + std::sqrt(a + 0.86595709258347853 * b + + 0.0048207809894260144 * c + 1.8611363115940525 * d + + 0.13886368840594743 * e + 0.12922212642709538 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.86595709258347853 * b + + 0.10890625570683385 * c + 1.8611363115940525 * d + + 0.66001895641514374 * e + 0.61419262306231814 * f) + + 0.056712962962962937 * + std::sqrt(a + 0.86595709258347853 * b + + 0.44888729929169013 * c + 1.8611363115940525 * d + + 1.3399810435848563 * e + 1.2469436885317342 * f) + + 0.03025074832140047 * + std::sqrt(a + 0.86595709258347853 * b + + 0.86595709258347853 * c + 1.8611363115940525 * d + + 1.8611363115940525 * e + 1.7319141851669571 * f); + } + + // the face is planar. then its area is 1/2 of the norm of the + // cross product of the two diagonals + const Tensor<1, 3> v12 = accessor.vertex(2) - accessor.vertex(1); + const Tensor<1, 3> twice_area = cross_product_3d(v03, v12); + return 0.5 * twice_area.norm(); } diff --git a/tests/codim_one/measure_of_distorted_codim_one_face.cc b/tests/codim_one/measure_of_distorted_codim_one_face.cc index efb6d30c99..342a64f100 100644 --- a/tests/codim_one/measure_of_distorted_codim_one_face.cc +++ b/tests/codim_one/measure_of_distorted_codim_one_face.cc @@ -43,10 +43,14 @@ test() std::vector> cells; SubCellData subcelldata; - double tol = 1e-15; - vertices.push_back(Point<3>{10, 56, 0}); - vertices.push_back(Point<3>{22, 1, 0}); - vertices.push_back(Point<3>{15, 44, 0}); + double tol = 1e-12; + // vertices.push_back(Point<3>{10, 56, 0}); + // vertices.push_back(Point<3>{22, 1, 0}); + // vertices.push_back(Point<3>{15, 44, 0}); + // vertices.push_back(Point<3>{1, 1, 1}); + vertices.push_back(Point<3>{0, 0, 1}); + vertices.push_back(Point<3>{1, 0, -10}); + vertices.push_back(Point<3>{0, 1, -1}); vertices.push_back(Point<3>{1, 1, 1}); cells.resize(1); @@ -62,7 +66,7 @@ test() DoFHandler<2, 3> dof_handler(tria); dof_handler.distribute_dofs(fe); - QGauss<2> quadrature_formula(1); + QGauss<2> quadrature_formula(4); MappingQ1<2, 3> mapping; FEValues<2, 3> fe_values(mapping, fe, quadrature_formula, update_JxW_values); @@ -90,7 +94,7 @@ int main() { initlog(); - deallog << std::setprecision(5); + deallog << std::setprecision(8); test(); -- 2.39.5