]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fixed bug in SphericalManifold for 3D case. 3175/head
authorLuca Heltai <luca.heltai@sissa.it>
Sat, 1 Oct 2016 16:42:30 +0000 (18:42 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Sun, 2 Oct 2016 10:20:28 +0000 (12:20 +0200)
When R>>1 and C != Point<spacedim>(), SphericalManifold
computed directions by subtracting the center twice.

doc/news/changes.h
source/grid/manifold_lib.cc
tests/manifold/spherical_manifold_04.cc [new file with mode: 0644]
tests/manifold/spherical_manifold_04.output [new file with mode: 0644]

index cfe7aa904d24169e4ba2d1fffe5e391a4f540a61..85e76bcf5ab7953b6ae2390246f512838bf4c30f 100644 (file)
@@ -384,15 +384,21 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+ <li> Fixed: SphericalManifold now behaves correctly also when R>>1
+ and the center is not the origin. 
+ <br>
+ (Luca Heltai, 2016/10/01)
+ </li>
+
  <li> Improved: Some parts of mesh refinement are now parallelized.
- <br/>
+ <br>
  (Wolfgang Bangerth, 2016/09/27)
  </li>
 
  <li> 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.
- <br/>
+ <br>
  (Jonathan Robey, 2016/09/21)
  </li>
 
index 1feaf897054e5a85fd042f22d271aac0b42fe19b..23bc0222698f1f282279eb3bcda55c6de935ba79 100644 (file)
@@ -276,14 +276,16 @@ get_new_point (const std::vector<Point<spacedim> > &vertices,
                const std::vector<double> &weights) const
 {
   double rho = 0.0;
-  Point<spacedim> candidate;
+  Tensor<1,spacedim> candidate;
   for (unsigned int i = 0; i<vertices.size(); i++)
     {
       rho += (vertices[i]-center).norm()*weights[i];
       candidate += (vertices[i]-center)*weights[i];
     }
+  // Unit norm direction.
+  candidate /= candidate.norm();
 
-  return center+(rho/(candidate-center).norm())*(candidate-center);
+  return center+rho*candidate;
 }
 
 // ============================================================
diff --git a/tests/manifold/spherical_manifold_04.cc b/tests/manifold/spherical_manifold_04.cc
new file mode 100644 (file)
index 0000000..e41f010
--- /dev/null
@@ -0,0 +1,102 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+// test Volume of a Ball
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/numerics/vector_tools.h>
+#include <deal.II/numerics/data_out.h>
+
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_in.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/fe_evaluation.h>
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+
+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<dim> center;
+  for (unsigned int d=0; d < dim; d++)
+    center[d] = d;
+
+  Triangulation<dim> triangulation;
+  DoFHandler<dim> dof_handler(triangulation);
+  FE_Q<dim> fe(fe_degree);
+  QGauss<dim> quadrature_formula(n_q_points_1d);
+
+  GridGenerator::hyper_ball (triangulation,
+                             center,
+                             R);
+  triangulation.set_all_manifold_ids_on_boundary(0);
+  static SphericalManifold<dim> surface_description(center);
+  triangulation.set_manifold (0, surface_description);
+  triangulation.refine_global(global_mesh_refinement_steps);
+  dof_handler.distribute_dofs (fe);
+  MappingQ<dim> mapping(fe_degree);
+
+  FEValues<dim> fe_values (mapping, fe, quadrature_formula,
+                           update_quadrature_points |
+                           update_JxW_values);
+
+  typename DoFHandler<dim>::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<n_q_points; ++q)
+        volume += fe_values.JxW (q);
+    }
+
+  deallog << "Volume:       " << volume  << std::endl
+          << "Exact volume: " << 4.0*numbers::PI *std::pow(R,3.0)/3. << std::endl;
+
+  dof_handler.clear ();
+}
+
+
+using namespace dealii;
+
+int main (int argc, char *argv[])
+{
+  initlog();
+  test(15);
+  return 0;
+}
diff --git a/tests/manifold/spherical_manifold_04.output b/tests/manifold/spherical_manifold_04.output
new file mode 100644 (file)
index 0000000..87df8eb
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::Volume:       14137.2
+DEAL::Exact volume: 14137.2

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.