From 5cf251dba4060600ebfb90f23b0b18668e98f163 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sat, 1 Oct 2016 18:42:30 +0200 Subject: [PATCH] Fixed bug in SphericalManifold for 3D case. When R>>1 and C != Point(), SphericalManifold computed directions by subtracting the center twice. --- doc/news/changes.h | 10 +- source/grid/manifold_lib.cc | 6 +- tests/manifold/spherical_manifold_04.cc | 102 ++++++++++++++++++++ tests/manifold/spherical_manifold_04.output | 3 + 4 files changed, 117 insertions(+), 4 deletions(-) create mode 100644 tests/manifold/spherical_manifold_04.cc create mode 100644 tests/manifold/spherical_manifold_04.output diff --git a/doc/news/changes.h b/doc/news/changes.h index cfe7aa904d..85e76bcf5a 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -384,15 +384,21 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: SphericalManifold now behaves correctly also when R>>1 + and the center is not the origin. +
    + (Luca Heltai, 2016/10/01) +
  2. +
  3. Improved: Some parts of mesh refinement are now parallelized. -
    +
    (Wolfgang Bangerth, 2016/09/27)
  4. Improved: MGSmootherBlock is now able to use the shared memory pool for temporary vector allocation. The constructor requiring an external memory allocation has therefore been deprecated. -
    +
    (Jonathan Robey, 2016/09/21)
  5. diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc index 1feaf89705..23bc022269 100644 --- a/source/grid/manifold_lib.cc +++ b/source/grid/manifold_lib.cc @@ -276,14 +276,16 @@ get_new_point (const std::vector > &vertices, const std::vector &weights) const { double rho = 0.0; - Point candidate; + Tensor<1,spacedim> candidate; for (unsigned int i = 0; i + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +using namespace dealii; + +void test (const double R) +{ + const unsigned int dim = 3; + const unsigned int global_mesh_refinement_steps = 4; + const unsigned int fe_degree = 2; + const unsigned int n_q_points_1d = 3; + + // derived + Point center; + for (unsigned int d=0; d < dim; d++) + center[d] = d; + + Triangulation triangulation; + DoFHandler dof_handler(triangulation); + FE_Q fe(fe_degree); + QGauss quadrature_formula(n_q_points_1d); + + GridGenerator::hyper_ball (triangulation, + center, + R); + triangulation.set_all_manifold_ids_on_boundary(0); + static SphericalManifold surface_description(center); + triangulation.set_manifold (0, surface_description); + triangulation.refine_global(global_mesh_refinement_steps); + dof_handler.distribute_dofs (fe); + MappingQ mapping(fe_degree); + + FEValues fe_values (mapping, fe, quadrature_formula, + update_quadrature_points | + update_JxW_values); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + const unsigned int n_q_points = quadrature_formula.size(); + + double volume = 0.; + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + for (unsigned int q=0; q