From 703d2e913d0dfd256bf327807f83358b9795c4f3 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 11 Dec 2020 21:14:05 +0100 Subject: [PATCH] Fix FlatManifold::normal_vector for triangular faces --- include/deal.II/grid/reference_cell.h | 65 +++++++++++++++++++++++++++ source/grid/manifold.cc | 27 +++++++---- 2 files changed, 84 insertions(+), 8 deletions(-) diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 9b99f6ef76..28a314d2b4 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -118,6 +118,71 @@ namespace ReferenceCell return table[dim][n_vertices]; } + + /** + * Compute the value of the $i$-th linear shape function at location $\xi$ for + * a given reference-cell type. + */ + template + inline double + d_linear_shape_function(const Type & reference_cell, + const Point & xi, + const unsigned int i) + { + if (reference_cell == get_hypercube(dim)) + return GeometryInfo::d_linear_shape_function(xi, i); + + if (reference_cell == + Type::Tri) // see also Simplex::ScalarPolynomial::compute_value + { + switch (i) + { + case 0: + return 1.0 - xi[0] - xi[1]; + case 1: + return xi[0]; + case 2: + return xi[1]; + } + } + + Assert(false, ExcNotImplemented()); + + return 0.0; + } + + /** + * Compute the gradient of the $i$-th linear shape function at location $\xi$ + * for a given reference-cell type. + */ + template + inline Tensor<1, dim> + d_linear_shape_function_gradient(const Type & reference_cell, + const Point & xi, + const unsigned int i) + { + if (reference_cell == get_hypercube(dim)) + return GeometryInfo::d_linear_shape_function_gradient(xi, i); + + if (reference_cell == + Type::Tri) // see also Simplex::ScalarPolynomial::compute_grad + { + switch (i) + { + case 0: + return Point(-1.0, -1.0); + case 1: + return Point(+1.0, +0.0); + case 2: + return Point(+0.0, +1.0); + } + } + + Assert(false, ExcNotImplemented()); + + return Point(+0.0, +0.0, +0.0); + } + namespace internal { /** diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index b839f903fd..b834e3f686 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -892,8 +892,19 @@ FlatManifold::normal_vector( const unsigned int facedim = dim - 1; Point xi; - for (unsigned int i = 0; i < facedim; ++i) - xi[i] = 1. / 2; + + const auto face_reference_cell_type = face->reference_cell_type(); + + if (face_reference_cell_type == ReferenceCell::get_hypercube(facedim)) + { + for (unsigned int i = 0; i < facedim; ++i) + xi[i] = 1. / 2; + } + else + { + for (unsigned int i = 0; i < facedim; ++i) + xi[i] = 1. / 3; + } const double eps = 1e-12; Tensor<1, spacedim> grad_F[facedim]; @@ -901,17 +912,17 @@ FlatManifold::normal_vector( while (true) { Point F; - for (const unsigned int v : GeometryInfo::vertex_indices()) - F += face->vertex(v) * - GeometryInfo::d_linear_shape_function(xi, v); + for (const unsigned int v : face->vertex_indices()) + F += face->vertex(v) * ReferenceCell::d_linear_shape_function( + face_reference_cell_type, xi, v); for (unsigned int i = 0; i < facedim; ++i) { grad_F[i] = 0; - for (const unsigned int v : GeometryInfo::vertex_indices()) + for (const unsigned int v : face->vertex_indices()) grad_F[i] += - face->vertex(v) * - GeometryInfo::d_linear_shape_function_gradient(xi, v)[i]; + face->vertex(v) * ReferenceCell::d_linear_shape_function_gradient( + face_reference_cell_type, xi, v)[i]; } Tensor<1, facedim> J; -- 2.39.5