From ccb56e7fb59f4a4892f27fcff876846fa631b213 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 7 May 2016 15:13:35 -0400 Subject: [PATCH] Fix a bug in the parallelepiped mesh generator. This function generated cells with negative volume due to an orientation issue. In particular, the function used to generate cells (see tests/grid/grid_parallelepiped.output) with vertices (given here in lexical order) 0 0 0 0.5 0.5 0.5 0.5 0 which results in a twisted cell, instead of 0 0 0.5 0 0.5 0.5 0 0.5 which is correct. --- doc/news/changes.h | 7 + source/grid/grid_generator.cc | 73 ++++---- tests/grid/grid_parallelepiped.output | 4 +- tests/grid/grid_parallelepiped_03.output | 204 +++++++++++------------ tests/grid/grid_parallelepiped_05.output | 16 +- 5 files changed, 164 insertions(+), 140 deletions(-) diff --git a/doc/news/changes.h b/doc/news/changes.h index 7718a459df..c62ea5a10d 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -249,6 +249,13 @@ inconvenience this causes. (Joscha Gedicke, 2016/05/07) +
  • Fixed: The function GridGenerator::subdivided_parallelepiped and its + variants could generate two dimensional meshes with cells that had negative + Jacobians. This has been fixed. +
    + (David Wells, 2016/05/07) +
  • +
  • New: Added function GridOut::write_mesh_per_processor_as_vtu. This allows the visualization of a parallel finite element mesh that can be separated into each processor's owned and ghost cells. It also allows for the visualization of each level diff --git a/source/grid/grid_generator.cc b/source/grid/grid_generator.cc index d1aafe9c33..341b82f339 100644 --- a/source/grid/grid_generator.cc +++ b/source/grid/grid_generator.cc @@ -856,26 +856,43 @@ namespace GridGenerator const std::vector &subdivisions, const bool colorize) { - if (subdivisions.size()==0) + // The provided edges and subdivisions may not be valid, so possibly + // recompute both of them. + std_cxx11::array, dim> compute_edges = edges; + std::vector compute_subdivisions = subdivisions; + if (compute_subdivisions.size() == 0) { - std::vector new_subdivisions(dim, 1); - subdivided_parallelepiped(tria, origin, edges, new_subdivisions, colorize); - return; + compute_subdivisions.resize(dim, 1); } - Assert(subdivisions.size()==dim, ExcMessage("")); - + Assert(compute_subdivisions.size()==dim, + ExcMessage("One subdivision must be provided for each dimension.")); // check subdivisions for (unsigned int i=0; i0, ExcInvalidRepetitions(subdivisions[i])); - Assert (edges[i].norm()>0, ExcMessage("Edges in subdivided_parallelepiped() must not be degenerate.")); + Assert (compute_subdivisions[i]>0, ExcInvalidRepetitions(subdivisions[i])); + Assert (edges[i].norm()>0, + ExcMessage("Edges in subdivided_parallelepiped() must not be degenerate.")); + } + + /* + * Verify that the vectors are oriented in a counter clockwise direction in 2D. + */ + if (dim == 2) + { + const double plane_normal = edges[0][0]*edges[1][1] - edges[0][1]*edges[1][0]; + if (plane_normal < 0.0) + { + std::swap(compute_edges[0], compute_edges[1]); + std::swap(compute_subdivisions[0], compute_subdivisions[1]); + } } + // Check corners do not overlap (unique) for (unsigned int i=0; i > cells (n_cells); // Create fixed ordering of switch (dim) { case 1: - for (unsigned int x=0; x