From bb6ce7d34bcbbf20b1099fa677ca858f65eaaa68 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 31 Jan 2017 18:59:29 +0100 Subject: [PATCH] Fix support points on quad for MappingQGeneric and StraightBoundary. --- source/fe/mapping_q_generic.cc | 54 +++-- tests/mappings/mapping_q_mixed_manifolds.cc | 201 ++++++++++++++++++ .../mappings/mapping_q_mixed_manifolds.output | 13 ++ 3 files changed, 254 insertions(+), 14 deletions(-) create mode 100644 tests/mappings/mapping_q_mixed_manifolds.cc create mode 100644 tests/mappings/mapping_q_mixed_manifolds.output diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index a7a48963d6..b9c6c65dc5 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -3775,9 +3775,9 @@ add_quad_support_points(const Triangulation<3,3>::cell_iterator &cell, { const unsigned int faces_per_cell = GeometryInfo<3>::faces_per_cell; - static const StraightBoundary<3> straight_boundary; // used if face quad at boundary or entirely in the interior of the domain std::vector > quad_points ((polynomial_degree-1)*(polynomial_degree-1)); + std::vector > tmp_points; // loop over all faces and collect points on them for (unsigned int face_no=0; face_no::cell_iterator &cell, ExcInternalError()); #endif - // ask the boundary/manifold object to return intermediate points on it - get_intermediate_points_on_object(face->get_manifold(), line_support_points, - face, quad_points); - - // in 3D, the orientation, flip and rotation of the face might not - // match what we expect here, namely the standard orientation. thus - // reorder points accordingly. since a Mapping uses the same shape - // function as an FE_Q, we can ask a FE_Q to do the reordering for us. - for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, - face_orientation, - face_flip, - face_rotation)]); + // On a quad, we have to check whether the manifold should determine the + // point distribution or rather a weighted sum should be created. This + // is the same logic as in the compute_mapping_support_points + // function below. + if (dynamic_cast *>(&face->get_manifold()) == NULL) + { + // ask the boundary/manifold object to return intermediate points on it + get_intermediate_points_on_object(face->get_manifold(), line_support_points, + face, quad_points); + + // in 3D, the orientation, flip and rotation of the face might not + // match what we expect here, namely the standard orientation. thus + // reorder points accordingly. since a Mapping uses the same shape + // function as an FE_Q, we can ask a FE_Q to do the reordering for us. + for (unsigned int i=0; iadjust_quad_dof_index_for_face_orientation(i, + face_orientation, + face_flip, + face_rotation)]); + } + else + { + // need to extract the points surrounding a quad from the points + // already computed. First get the 4 vertices and then the points on + // the four lines + tmp_points.resize(4 + 4*(polynomial_degree-1)); + for (unsigned int v=0; v::vertices_per_cell; ++v) + tmp_points[v] = a[GeometryInfo<3>::face_to_cell_vertices(face_no,v)]; + if (polynomial_degree > 1) + for (unsigned int line=0; line::lines_per_cell; ++line) + for (unsigned int i=0; i::vertices_per_cell + + (polynomial_degree-1)* + GeometryInfo<3>::face_to_cell_lines(face_no,line) + i]; + add_weighted_interior_points (support_point_weights_on_quad, tmp_points); + a.insert(a.end(), tmp_points.begin()+4+4*(polynomial_degree-1), + tmp_points.end()); + } } } diff --git a/tests/mappings/mapping_q_mixed_manifolds.cc b/tests/mappings/mapping_q_mixed_manifolds.cc new file mode 100644 index 0000000000..fa0840c718 --- /dev/null +++ b/tests/mappings/mapping_q_mixed_manifolds.cc @@ -0,0 +1,201 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Shows the shape functions implemented and computes the area of cells. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +const double D = 0.1; +const double R = D/2.0; +const double R_1 = 1.2*R; +const double R_2 = 1.7*R; +const double H = 0.41; +const double X_0 = 0.0; +const double X_1 = 0.3; +const double X_C = 0.5; // center +const double X_2 = 0.7; + +const double Y_0 = 0.0; +const double Y_C = 0.2; // center + +const unsigned int MANIFOLD_ID = 1; + + +void create_triangulation(Triangulation<2> &tria) +{ + AssertThrow(std::abs((X_2-X_1) - 2.0*(X_C-X_1))<1.0e-12, ExcMessage("Geometry parameters X_1,X_2,X_C invalid!")); + SphericalManifold<2> spherical_manifold(Point<2>(X_C,Y_C)); + + Triangulation<2> circle_1, circle_2, circle_tmp, middle, middle_tmp, middle_tmp2, tmp_3D; + std::vector ref_1(2, 2); + ref_1[1] = 2; + + // create middle part first as a hyper shell + const double outer_radius = (X_2-X_1)/2.0; + const unsigned int n_cells = 4; + GridGenerator::hyper_shell(middle, Point<2>(X_C, Y_C), R_2, outer_radius, n_cells, true); + middle.set_all_manifold_ids(MANIFOLD_ID); + middle.set_manifold(MANIFOLD_ID, spherical_manifold); + middle.refine_global(1); + + // two inner circles in order to refine towards the cylinder surface + const unsigned int n_cells_circle = 8; + GridGenerator::hyper_shell(circle_1, Point<2>(X_C, Y_C), R, R_1, n_cells_circle, true); + circle_1.set_all_manifold_ids(MANIFOLD_ID); + circle_1.set_manifold(MANIFOLD_ID,spherical_manifold); + + GridGenerator::hyper_shell(circle_2, Point<2>(X_C, Y_C), R_1, R_2, n_cells_circle, true); + circle_2.set_all_manifold_ids(MANIFOLD_ID); + circle_2.set_manifold(MANIFOLD_ID,spherical_manifold); + + // then move the vertices to the points where we want them to be to create a slightly asymmetric cube with a hole + for (Triangulation<2>::cell_iterator cell = middle.begin(); cell != middle.end(); ++cell) + { + for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v) + { + Point<2> &vertex = cell->vertex(v); + if (std::abs(vertex[0] - X_2) < 1e-10 && std::abs(vertex[1] - Y_C) < 1e-10) + { + vertex = Point<2>(X_2, H/2.0); + } + else if (std::abs(vertex[0] - (X_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) + { + vertex = Point<2>(X_2, H); + } + else if (std::abs(vertex[0] - (X_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) + { + vertex = Point<2>(X_2, Y_0); + } + else if (std::abs(vertex[0] - X_C) < 1e-10 && std::abs(vertex[1] - (Y_C +(X_2-X_1)/2.0)) < 1e-10) + { + vertex = Point<2>(X_C, H); + } + else if (std::abs(vertex[0] - X_C) < 1e-10 && std::abs(vertex[1] - (Y_C-(X_2-X_1)/2.0)) < 1e-10) + { + vertex = Point<2>(X_C, Y_0); + } + else if (std::abs(vertex[0] - (X_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C + (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) + { + vertex = Point<2>(X_1, H); + } + else if (std::abs(vertex[0] - (X_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10 && std::abs(vertex[1] - (Y_C - (X_2-X_1)/2.0/std::sqrt(2))) < 1e-10) + { + vertex = Point<2>(X_1, Y_0); + } + else if (std::abs(vertex[0] - X_1) < 1e-10 && std::abs(vertex[1] - Y_C) < 1e-10) + { + vertex = Point<2>(X_1, H/2.0); + } + } + } + + // must copy the triangulation because we cannot merge triangulations with refinement... + GridGenerator::flatten_triangulation(middle, middle_tmp); + GridGenerator::merge_triangulations(circle_1,circle_2,circle_tmp); + GridGenerator::merge_triangulations(middle_tmp,circle_tmp,tria); + + // Set the cylinder boundary to 2, outflow to 1, the rest to 0. + //tria.set_all_manifold_ids(0); + for (Triangulation<2>::active_cell_iterator cell=tria.begin(); cell != tria.end(); ++cell) + { + if (Point<2>(X_C,Y_C).distance(cell->center())<= R_2) + cell->set_all_manifold_ids(MANIFOLD_ID); + } +} + +void create_triangulation(Triangulation<3> &tria) +{ + Triangulation<2> tria_2d; + create_triangulation(tria_2d); + GridGenerator::extrude_triangulation(tria_2d, 3, H, tria); + + // Set the cylinder boundary to 2, outflow to 1, the rest to 0. + tria.set_all_manifold_ids(0); + for (Triangulation<3>::active_cell_iterator cell=tria.begin(); cell != tria.end(); ++cell) + { + if (Point<3>(X_C,Y_C,cell->center()[2]).distance(cell->center())<= R_2) + cell->set_all_manifold_ids(MANIFOLD_ID); + } +} + + + +template +void test() +{ + Point center; + center[0] = X_C; + center[1] = Y_C; + Point direction; + direction[dim-1] = 1.; + + std_cxx11::shared_ptr > cylinder_manifold + (dim == 2 ? static_cast*>(new SphericalManifold(center)) : + static_cast*>(new CylindricalManifold(direction, center))); + Triangulation tria; + create_triangulation(tria); + tria.set_manifold(MANIFOLD_ID, *cylinder_manifold); + + FE_Nothing fe; + for (unsigned int degree = 1; degree < 7; ++degree) + { + MappingQGeneric mapping(degree); + QGauss quad(degree+1); + FEValues fe_values(mapping, fe, quad, update_JxW_values); + double sum = 0.; + for (typename Triangulation::active_cell_iterator cell=tria.begin_active(); + cell != tria.end(); ++cell) + { + fe_values.reinit(cell); + double local_sum = 0; + for (unsigned int q=0; q(); + test<3>(); +} diff --git a/tests/mappings/mapping_q_mixed_manifolds.output b/tests/mappings/mapping_q_mixed_manifolds.output new file mode 100644 index 0000000000..19c6438ed5 --- /dev/null +++ b/tests/mappings/mapping_q_mixed_manifolds.output @@ -0,0 +1,13 @@ + +DEAL::Volume 2D mapping degree 1: 0.156929 error: 0.00501399 +DEAL::Volume 2D mapping degree 2: 0.156152 error: 3.91474e-05 +DEAL::Volume 2D mapping degree 3: 0.156147 error: 6.04868e-06 +DEAL::Volume 2D mapping degree 4: 0.156146 error: 5.88206e-08 +DEAL::Volume 2D mapping degree 5: 0.156146 error: 8.80543e-09 +DEAL::Volume 2D mapping degree 6: 0.156146 error: 8.99033e-11 +DEAL::Volume 3D mapping degree 1: 0.0643409 error: 0.00501399 +DEAL::Volume 3D mapping degree 2: 0.0640224 error: 3.91474e-05 +DEAL::Volume 3D mapping degree 3: 0.0640203 error: 6.04868e-06 +DEAL::Volume 3D mapping degree 4: 0.0640199 error: 5.88206e-08 +DEAL::Volume 3D mapping degree 5: 0.0640199 error: 8.80543e-09 +DEAL::Volume 3D mapping degree 6: 0.0640199 error: 8.99045e-11 -- 2.39.5