* cells; again this class is then equivalent to using MappingQGeneric right
* away.
*
+ * <h4>Behavior along curved boundaries and with different manifolds</h4>
+ *
+ * As described above, one often only knows a manifold description of a
+ * surface but not the interior of the computational domain. Thus, a
+ * StraightBoundary object will be assigned to the interior entities that
+ * describes a usual planar coordinate system where the additional points of
+ * this class are placed exactly according to a bi-/trilinear mapping. When
+ * combined with a non-flat manifold on the boundary, for example a circle,
+ * the two manifold descriptions are in general incompatible. For example, a
+ * StraightBoundary would put an interior point located at some small distance
+ * epsilon away from the boundary along a flat line and thus in general
+ * outside the concave part of a circle. If the polynomial degree of
+ * MappingQGeneric is sufficiently high, the transformation from the reference
+ * cell to such a cell would in general contain inverted regions close to the
+ * boundary.
+ *
+ * In order to avoid this situation, this class applies a smoothing on cells
+ * adjacent to the boundary by using so-called Laplace smoothing. In the
+ * algorithm computing new points, all the entities of the cells are passed
+ * through hierarchically, starting from the lines to the quads and finally
+ * hexes. The elements higher up in the hierarchy that sit on a
+ * StraightBoundary will then get their points interpolated from all the
+ * surrounding points and not just the corner points. If only a line is
+ * associated to a curved boundary but the adjacent quad is on a flat
+ * manifold, the points inside the quad will be computed according to the
+ * deformed line and thus always result in a well-defined transformation. If
+ * this smoothing is undesired, the optional argument @p
+ * smooth_support_points_on_flat_manifold can be used to disable the
+ * smoothing, placing the additional points strictly according to the
+ * manifold.
+ *
+ * While this smoothing approach works well for low and medium convergence
+ * orders up to approximately three to four, this mechanism has an inherent
+ * shortcoming because it switches from a curved manifold to a flat manifold
+ * within one layer of elements. This will cause the Jacobian transformation
+ * to have jumps between the first and second element layer that can reduce
+ * the order of convergence. For example, the convergence rates for solving
+ * the Laplacian on a circle where only the boundary is deformed and the above
+ * mesh smoothing algorithm is applied will typically not exceed 3.5, even for
+ * fourth or fifth degree polynomials. In such a case, the curved manifold
+ * needs to be switched to a flat manifold in a smooth way that does not
+ * depend on the mesh size and eventuell covers a whole layer of cells.
+ *
* @author Ralf Hartmann, 2000, 2001, 2005; Guido Kanschat 2000, 2001,
* Wolfgang Bangerth, 2015
*/
* The value of @p use_mapping_q_on_all_cells is ignored if @p dim is not
* equal to @p spacedim, i.e., if we are considering meshes on surfaces
* embedded into higher dimensional spaces.
+ *
+ * The optional parameter @p smooth_support_points_on_flat_manifold controls
+ * whether smoothing adjancent to StraightBoundary objects according to the
+ * general class description should be enabled (default) or not.
*/
MappingQ (const unsigned int polynomial_degree,
- const bool use_mapping_q_on_all_cells = false);
+ const bool use_mapping_q_on_all_cells = false,
+ const bool smooth_support_points_on_flat_manifold = true);
/**
* Copy constructor.
* a manifold object to all cells that can then also be used by the mapping
* classes for higher order mappings.
*
+ * <h4>Behavior along curved boundaries and with different manifolds</h4>
+ *
+ * As described above, one often only knows a manifold description of a
+ * surface but not the interior of the computational domain. Thus, a
+ * StraightBoundary object will be assigned to the interior entities that
+ * describes a usual planar coordinate system where the additional points of
+ * this class are placed exactly according to a bi-/trilinear mapping. When
+ * combined with a non-flat manifold on the boundary, for example a circle,
+ * the two manifold descriptions are in general incompatible. For example, a
+ * StraightBoundary would put an interior point located at some small distance
+ * epsilon away from the boundary along a flat line and thus in general
+ * outside the concave part of a circle. If the polynomial degree of
+ * MappingQGeneric is sufficiently high, the transformation from the reference
+ * cell to such a cell would in general contain inverted regions close to the
+ * boundary.
+ *
+ * In order to avoid this situation, this class applies a smoothing on cells
+ * adjacent to the boundary by using so-called Laplace smoothing. In the
+ * algorithm computing new points, all the entities of the cells are passed
+ * through hierarchically, starting from the lines to the quads and finally
+ * hexes. The elements higher up in the hierarchy that sit on a
+ * StraightBoundary will then get their points interpolated from all the
+ * surrounding points and not just the corner points. If only a line is
+ * associated to a curved boundary but the adjacent quad is on a flat
+ * manifold, the points inside the quad will be computed according to the
+ * deformed line and thus always result in a well-defined transformation. If
+ * this smoothing is undesired, the optional argument @p
+ * smooth_support_points_on_flat_manifold can be used to disable the
+ * smoothing, placing the additional points strictly according to the
+ * manifold.
+ *
+ * While this smoothing approach works well for low and medium convergence
+ * orders up to approximately three to four, this mechanism has an inherent
+ * shortcoming because it switches from a curved manifold to a flat manifold
+ * within one layer of elements. This will cause the Jacobian transformation
+ * to have jumps between the first and second element layer that can reduce
+ * the order of convergence. For example, the convergence rates for solving
+ * the Laplacian on a circle where only the boundary is deformed and the above
+ * mesh smoothing algorithm is applied will typically not exceed 3.5, even for
+ * fourth or fifth degree polynomials. In such a case, the curved manifold
+ * needs to be switched to a flat manifold in a smooth way that does not
+ * depend on the mesh size and eventuell covers a whole layer of cells.
*
* @author Wolfgang Bangerth, 2015
*/
* Constructor. @p polynomial_degree denotes the polynomial degree of the
* polynomials that are used to map cells from the reference to the real
* cell.
+ *
+ * The optional parameter @p smooth_support_points_on_flat_manifold controls
+ * whether smoothing adjancent to StraightBoundary objects according to the
+ * general class description should be enabled (default) or not.
*/
- MappingQGeneric (const unsigned int polynomial_degree);
+ MappingQGeneric (const unsigned int polynomial_degree,
+ const bool smooth_support_points_on_flat_manifold = true);
/**
* Copy constructor.
*/
const unsigned int polynomial_degree;
+ /**
+ * Stores whether we want to smooth the placement of interior points on flat
+ * manifolds.
+ */
+ const bool smooth_support_points_on_flat_manifold;
+
/*
* The default line support points. These are used when computing
* the location in real space of the support points on lines and
template<int dim, int spacedim>
MappingQ<dim,spacedim>::MappingQ (const unsigned int degree,
- const bool use_mapping_q_on_all_cells)
+ const bool use_mapping_q_on_all_cells,
+ const bool smooth_support_points_on_flat_manifold)
:
polynomial_degree (degree),
// created via the shared_ptr objects
qp_mapping (this->polynomial_degree>1
?
- std_cxx11::shared_ptr<const MappingQGeneric<dim,spacedim> >(new MappingQGeneric<dim,spacedim>(degree))
+ std_cxx11::shared_ptr<const MappingQGeneric<dim,spacedim> >(new MappingQGeneric<dim,spacedim>(degree, smooth_support_points_on_flat_manifold))
:
q1_mapping)
{}
template<int dim, int spacedim>
-MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p)
+MappingQGeneric<dim,spacedim>::MappingQGeneric (const unsigned int p,
+ const bool smooth_support_points_on_flat_manifold)
:
polynomial_degree(p),
+ smooth_support_points_on_flat_manifold(smooth_support_points_on_flat_manifold),
line_support_points(this->polynomial_degree+1),
fe_q(dim == 3 ? new FE_Q<dim>(this->polynomial_degree) : 0),
support_point_weights_on_quad (compute_support_point_weights_on_quad<dim>(this->polynomial_degree)),
MappingQGeneric<dim,spacedim>::MappingQGeneric (const MappingQGeneric<dim,spacedim> &mapping)
:
polynomial_degree(mapping.polynomial_degree),
+ smooth_support_points_on_flat_manifold(mapping.smooth_support_points_on_flat_manifold),
line_support_points(mapping.line_support_points),
fe_q(dim == 3 ? new FE_Q<dim>(*mapping.fe_q) : 0),
support_point_weights_on_quad (mapping.support_point_weights_on_quad),
// 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<const Boundary<3,3> *>(&face->get_manifold()) == NULL)
+ // is the same logic as in the compute_mapping_support_points function
+ // below: if we have a StraightBoundary but no "real" boundary
+ // description doing something else, we should go for the smoothing of
+ // the support points through add_weighted_interior_points
+ if (smooth_support_points_on_flat_manifold == false ||
+ std::string(typeid(face->get_manifold()).name()).find("StraightBoundary") ==
+ std::string::npos)
{
// ask the boundary/manifold object to return intermediate points on it
get_intermediate_points_on_object(face->get_manifold(), line_support_points,
// (outer) points
add_line_support_points(cell, a);
- // then get the support points on the quad if we are on a
- // manifold, otherwise compute them from the points around it
+ // then get the support points on the quad. if we have a
+ // StraightBoundary but no "real" boundary description doing something
+ // else, we should go for the smoothing of the support points through
+ // add_weighted_interior_points. we need to go through typeid because
+ // curved boundaries are derived from StraightBoundary and thus a
+ // plain dynamic_cast will not suffice.
if (dim != spacedim)
add_quad_support_points(cell, a);
- // TODO: use get_intermediate_points_on_object for the else case
- // unconditionally as soon as the interface of Boundary is fixed
- else if (dynamic_cast<const Boundary<dim,spacedim> *>(&cell->get_manifold()) == NULL)
+ else if (smooth_support_points_on_flat_manifold == false ||
+ std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") ==
+ std::string::npos)
{
std::vector<Point<spacedim> > quad_points (Utilities::fixed_power<dim>(polynomial_degree-1));
get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, quad_points);
add_line_support_points (cell, a);
add_quad_support_points (cell, a);
- // then compute the interior points
- // TODO: remove else case as soon as boundary is fixed
- if (dynamic_cast<const Boundary<dim,spacedim> *>(&cell->get_manifold()) == NULL)
+ // then compute the interior points. if we have a StraightBoundary but
+ // no "real" boundary description doing something else, we should go
+ // for the smoothing of the support points through
+ // add_weighted_interior_points. we need to go through typeid because
+ // curved boundaries are derived from StraightBoundary and thus a
+ // plain dynamic_cast will not suffice.
+ if (smooth_support_points_on_flat_manifold == false ||
+ std::string(typeid(cell->get_manifold()).name()).find("StraightBoundary") ==
+ std::string::npos)
{
std::vector<Point<spacedim> > hex_points (Utilities::fixed_power<dim>(polynomial_degree-1));
get_intermediate_points_on_object(cell->get_manifold(), line_support_points, cell, hex_points);
+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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 <deal.II/base/quadrature_lib.h>
-#include <deal.II/base/logstream.h>
-#include <deal.II/lac/vector.h>
-#include <deal.II/grid/tria.h>
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/grid/grid_generator.h>
-#include <deal.II/grid/grid_out.h>
-#include <deal.II/grid/tria_boundary_lib.h>
-#include <deal.II/fe/mapping_q_generic.h>
-#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe_nothing.h>
-#include <deal.II/grid/manifold_lib.h>
-#include <vector>
-#include <fstream>
-#include <iomanip>
-#include <string>
-#include <sstream>
-
-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<unsigned int> 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 <int dim>
-void test()
-{
- Point<dim> center;
- center[0] = X_C;
- center[1] = Y_C;
- Point<dim> direction;
- direction[dim-1] = 1.;
-
- std_cxx11::shared_ptr<Manifold<dim> > cylinder_manifold
- (dim == 2 ? static_cast<Manifold<dim>*>(new SphericalManifold<dim>(center)) :
- static_cast<Manifold<dim>*>(new CylindricalManifold<dim>(direction, center)));
- Triangulation<dim> tria;
- create_triangulation(tria);
- tria.set_manifold(MANIFOLD_ID, *cylinder_manifold);
-
- FE_Nothing<dim> fe;
- for (unsigned int degree = 1; degree < 7; ++degree)
- {
- MappingQGeneric<dim> mapping(degree);
- QGauss<dim> quad(degree+1);
- FEValues<dim> fe_values(mapping, fe, quad, update_JxW_values);
- double sum = 0.;
- for (typename Triangulation<dim>::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<quad.size(); ++q)
- local_sum += fe_values.JxW(q);
- sum += local_sum;
- }
- const double reference = (dim==2 ? 1. : H) * (H*(X_2-X_1)-D*D*numbers::PI*0.25);
- deallog << "Volume " << dim << "D mapping degree " << degree << ": "
- << sum << " error: " << (sum-reference)/reference << std::endl;
- }
-}
-
-
-int main()
-{
- std::ofstream logfile ("output");
- deallog.attach(logfile);
- deallog.threshold_double(1.e-14);
-
- test<2>();
- test<3>();
-}
projects / dealii.git / commitdiff