From: Nicola Giuliani Date: Wed, 8 Jan 2020 14:40:26 +0000 (+0100) Subject: added measure computation for distorted 2d cells X-Git-Tag: v9.2.0-rc1~160^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9260%2Fhead;p=dealii.git added measure computation for distorted 2d cells --- diff --git a/source/grid/tria_accessor.cc b/source/grid/tria_accessor.cc index 290917877c..62919885f6 100644 --- a/source/grid/tria_accessor.cc +++ b/source/grid/tria_accessor.cc @@ -1286,42 +1286,35 @@ namespace // a 2d face in 3d space + template double - measure(const dealii::TriaAccessor<2, 3, 3> &accessor) + measure(const dealii::TriaAccessor<2, dim, 3> &accessor) { - // 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) - { - Assert( - false, - ExcMessage( - "Computing the measure of a nonplanar face is not implemented!")); - return std::numeric_limits::quiet_NaN(); - } + // In general the area can be computed as + // 0.25*(v_0+v_1-v_2-v_3)*(v_0-v_1+v_2-v_3) - // 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(); + const Tensor<1, 3> piece_1 = accessor.vertex(0) + accessor.vertex(1) - + accessor.vertex(2) - accessor.vertex(3); + const Tensor<1, 3> piece_2 = accessor.vertex(0) - accessor.vertex(1) + + accessor.vertex(2) - accessor.vertex(3); + + return 0.25 * cross_product_3d(piece_1, piece_2).norm(); } + // // a 2d cell in 3d space + // double + // measure(const dealii::TriaAccessor<2, 2, 3> &accessor) + // { + // // In general the area can be computed as + // // 0.25*(v_0+v_1-v_2-v_3)*(v_0-v_1+v_2-v_3) + + // const Tensor<1, 3> piece_1 = accessor.vertex(0) + accessor.vertex(1) - + // accessor.vertex(2) - accessor.vertex(3); + // const Tensor<1, 3> piece_2 = accessor.vertex(0) - accessor.vertex(1) + + // accessor.vertex(2) - accessor.vertex(3); + // return 0.25 * cross_product_3d(piece_1, piece_2).norm(); + // } + template diff --git a/tests/codim_one/measure_of_distorted_codim_one_face.cc b/tests/codim_one/measure_of_distorted_codim_one_face.cc new file mode 100644 index 0000000000..efb6d30c99 --- /dev/null +++ b/tests/codim_one/measure_of_distorted_codim_one_face.cc @@ -0,0 +1,97 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2014 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test by Nicola Giuliani: compute the measure of a 2d distorted cell in codim +// one + +#include + +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "../tests.h" + + + +void +test() +{ + Triangulation<2, 3> tria; + std::vector> vertices; + 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}); + vertices.push_back(Point<3>{1, 1, 1}); + + cells.resize(1); + + cells[0].vertices[0] = 0; + cells[0].vertices[1] = 1; + cells[0].vertices[2] = 2; + cells[0].vertices[3] = 3; + + tria.create_triangulation(vertices, cells, subcelldata); + + FE_Q<2, 3> fe(1); + DoFHandler<2, 3> dof_handler(tria); + dof_handler.distribute_dofs(fe); + + QGauss<2> quadrature_formula(1); + MappingQ1<2, 3> mapping; + FEValues<2, 3> fe_values(mapping, fe, quadrature_formula, update_JxW_values); + + double sum_1 = 0; + double sum_2 = 0; + auto cell = dof_handler.begin_active(); + auto accessor = tria.begin_active(); + + sum_1 = cell->measure(); + fe_values.reinit(cell); + for (unsigned int q = 0; q < quadrature_formula.size(); ++q) + { + sum_2 += fe_values.JxW(q); + } + if (std::abs(sum_1 - sum_2) < tol) + deallog << "Test passed" << std::endl; + else + deallog << sum_1 << " " << sum_2 << std::endl; + deallog << std::endl; +} + + + +int +main() +{ + initlog(); + deallog << std::setprecision(5); + + + test(); +} diff --git a/tests/codim_one/measure_of_distorted_codim_one_face.output b/tests/codim_one/measure_of_distorted_codim_one_face.output new file mode 100644 index 0000000000..7ca18d1d8c --- /dev/null +++ b/tests/codim_one/measure_of_distorted_codim_one_face.output @@ -0,0 +1,3 @@ + +DEAL::Test passed +DEAL::