]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Added buggy spherical manifold.
authorLuca Heltai <luca.heltai@sissa.it>
Fri, 1 Aug 2014 09:08:04 +0000 (11:08 +0200)
committerLuca Heltai <luca.heltai@gmail.com>
Wed, 6 Aug 2014 09:20:15 +0000 (11:20 +0200)
Fixed compilation bug on Mac.

Modified tests for spherical manifold.

12 files changed:
include/deal.II/grid/manifold_lib.h [new file with mode: 0644]
source/grid/CMakeLists.txt
source/grid/manifold.cc
source/grid/manifold_lib.cc [new file with mode: 0644]
source/grid/manifold_lib.inst.in [new file with mode: 0644]
source/grid/tria_boundary_lib.cc
tests/manifold/spherical_manifold_01.cc [new file with mode: 0644]
tests/manifold/spherical_manifold_01.output [new file with mode: 0644]
tests/manifold/spherical_manifold_02.cc [new file with mode: 0644]
tests/manifold/spherical_manifold_02.output [new file with mode: 0644]
tests/manifold/spherical_manifold_03.cc [new file with mode: 0644]
tests/manifold/spherical_manifold_03.output [new file with mode: 0644]

diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h
new file mode 100644 (file)
index 0000000..6a1f7ee
--- /dev/null
@@ -0,0 +1,87 @@
+// ---------------------------------------------------------------------
+// $Id: tria_boundary_lib.h 30130 2013-07-23 13:01:18Z heltai $
+//
+// Copyright (C) 1999 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef __deal2__manifold_lib_h
+#define __deal2__manifold_lib_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/grid/manifold.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+/**
+ * Manifold description for a spherical space coordinate system. 
+ *
+ * You can use this Manifold object to describe any sphere, circle,
+ * hypersphere or hyperdisc in two or three dimensions, both as a
+ * co-dimension one manifold descriptor or as co-dimension zero
+ * manifold descriptor.
+ *
+ * The two template arguments match the meaning of the two template
+ * arguments in Triangulation<dim, spacedim>, however this Manifold
+ * can be used to describe both thin and thick objects, and the
+ * behavior is identical when dim <= spacedim, i.e., the functionality
+ * of SphericalManifold<2,3> is identical to SphericalManifold<3,3>.
+ *
+ * @ingroup manifold
+ *
+ * @author Luca Heltai, 2014
+ */
+template <int dim, int spacedim>
+class SphericalManifold : public ManifoldChart<dim, spacedim, spacedim> 
+{
+public:
+  /**
+   * The Constructor takes the center of the spherical coordinates
+   * system. This class uses the pull_back and push_forward mechanism
+   * to transform from cartesian to spherical coordinate systems,
+   * taking into account the periodicity of base Manifold.
+   */
+  SphericalManifold(const Point<spacedim> center = Point<spacedim>());
+
+  /**
+   * Pull back the given point from the Euclidean space. Will return
+   * the polar coordinates associated with the point @p space_point.
+   */
+  virtual Point<spacedim>
+  pull_back(const Point<spacedim> &space_point) const;
+
+  /**
+   * Given a point in the spherical coordinate system, this method
+   * returns the Euclidean coordinates associated to the polar
+   * coordinates @p chart_point.
+   */
+  virtual Point<spacedim>
+  push_forward(const Point<spacedim> &chart_point) const;
+  
+  
+  /**
+   * The center of the spherical coordinate system.
+   */
+  const Point<spacedim> center;
+private:
+  
+  /** Helper function which returns the periodicity associated with
+      this coordinate system, according to dim, chartdim, and
+      spacedim. */
+  static Point<spacedim> get_periodicity();
+};
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 1b074c95235cae4b3d51a6821b3c49004f4a683d..18e30fe8c7d52b9b46407b28aaa87908a9c32288 100644 (file)
@@ -25,6 +25,7 @@ SET(_src
   grid_tools.cc
   intergrid_map.cc
   manifold.cc
+  manifold_lib.cc
   persistent_tria.cc
   tria_accessor.cc
   tria_boundary.cc
index 1685bbc1bec7f5d8c3a1f8c606d40bb5c0526fa4..6d4dd06c81359fbab34909a71dd86dd8560d05fc 100644 (file)
@@ -64,7 +64,7 @@ namespace Manifolds {
   
   template <typename OBJECT, int spacedim>
   void get_default_quadrature(const OBJECT& obj, 
-                             Quadrature<spacedim> &quad, bool with_laplace = false
+                             Quadrature<spacedim> &quad, bool with_laplace) 
   {
     std::vector<Point<spacedim> > sp;
     std::vector<double> wp;
diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc
new file mode 100644 (file)
index 0000000..d410187
--- /dev/null
@@ -0,0 +1,113 @@
+// ---------------------------------------------------------------------
+// $Id: manifold_lib.cc 30130 2013-07-23 13:01:18Z heltai $
+//
+// Copyright (C) 2013 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.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/base/tensor.h>
+#include <cmath>
+
+DEAL_II_NAMESPACE_OPEN
+
+template <int dim, int spacedim>
+SphericalManifold<dim,spacedim>::SphericalManifold(const Point<spacedim> center):
+  ManifoldChart<dim,spacedim,spacedim>(SphericalManifold<dim,spacedim>::get_periodicity()),
+  center(center)
+{
+  Assert(spacedim != 1, ExcImpossibleInDim(1));
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+SphericalManifold<dim,spacedim>::get_periodicity()  {
+  Point<spacedim> periodicity;
+  periodicity[spacedim-1] = 2*numbers::PI; // theta and phi period.
+  return periodicity;
+}
+
+
+
+template <int dim, int spacedim>
+Point<spacedim>
+SphericalManifold<dim,spacedim>::push_forward(const Point<spacedim> &spherical_point) const {
+  Assert(spherical_point[0] >=0.0, 
+        ExcMessage("Negative radius for given point."));
+  const double &rho = spherical_point[0];
+  const double &theta = spherical_point[1];
+  
+  Point<spacedim> p;
+  if(rho > 1e-10) 
+    switch(spacedim) {
+    case 2: 
+      p[0] = rho*cos(theta);
+      p[1] = rho*sin(theta);
+      break;
+    case 3:
+      {
+       const double &phi= spherical_point[2];
+       p[0] = rho*sin(theta)*cos(phi);
+       p[1] = rho*sin(theta)*sin(phi);
+       p[2] = rho*cos(theta);
+      }
+      break;
+    default:
+      Assert(false, ExcInternalError());
+  }
+  return p+center;
+}
+
+template <int dim, int spacedim>
+Point<spacedim>
+SphericalManifold<dim,spacedim>::pull_back(const Point<spacedim> &space_point) const {
+  const Point<spacedim> R = space_point-center;
+  const double rho = R.norm();
+  const double &x = R[0];
+  const double &y = R[1];
+  
+  Point<spacedim> p;
+  p[0] = rho;
+  
+  switch(spacedim) {
+  case 2: 
+    p[1] = atan2(y,x);
+    if(p[1] < 0)
+      p[1] += 2*numbers::PI;
+    break;
+  case 3:
+    {
+      const double &z = R[2];
+      p[2] = atan2(y,x); // phi
+      if(p[2] < 0)
+       p[2] += 2*numbers::PI; // phi is periodic
+      p[1] = atan(sqrt(x*x+y*y)/z)+numbers::PI;  // theta
+    }
+    break;
+  default:
+    Assert(false, ExcInternalError());
+  }
+  return p;
+}
+               
+// explicit instantiations
+template class SphericalManifold<1,2>;
+template class SphericalManifold<2,2>;
+template class SphericalManifold<2,3>;
+template class SphericalManifold<3,3>;
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/grid/manifold_lib.inst.in b/source/grid/manifold_lib.inst.in
new file mode 100644 (file)
index 0000000..11bf419
--- /dev/null
@@ -0,0 +1,25 @@
+// ---------------------------------------------------------------------
+// $Id: tria_boundary_lib.inst.in 30130 2013-07-23 13:01:18Z heltai $
+//
+// Copyright (C) 1999 - 2013 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.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS)
+  {
+    template class PolarManifold<deal_II_dimension>;
+    template class CylinderManifold<deal_II_dimension>;
+//     template class ConeBoundary<deal_II_dimension>;
+  }
+
+
index 0043e63e3b87796f7a3e55910e945cc540732c92..a2ff2f43df6d1c4d2151b8c907bae46b7690e7be 100644 (file)
@@ -96,7 +96,7 @@ Point<3>
 CylinderBoundary<3>::
 get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
 {
-  const Point<3> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+  const Point<3> middle = StraightBoundary<3,3>::get_new_point_on_quad (quad);
 
   // same algorithm as above
   const unsigned int spacedim = 3;
@@ -355,7 +355,7 @@ get_new_point_on_quad (const Triangulation<3>::quad_iterator &quad) const
 
   const Point<dim> axis = x_1 - x_0;
   // Compute the middle point of the quad.
-  const Point<dim> middle = StraightBoundary<3>::get_new_point_on_quad (quad);
+  const Point<dim> middle = StraightBoundary<3,3>::get_new_point_on_quad (quad);
   // Same algorithm as above: To project it on the boundary of the cone we
   // first compute the orthogonal projection of the middle point onto the axis
   // of the cone.
diff --git a/tests/manifold/spherical_manifold_01.cc b/tests/manifold/spherical_manifold_01.cc
new file mode 100644 (file)
index 0000000..cf36d3b
--- /dev/null
@@ -0,0 +1,68 @@
+//----------------------------  spherical_manifold_01.cc  ---------------------------
+//    Copyright (C) 2011, 2013 by the mathLab team.
+//
+//    This file is subject to LGPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  spherical_manifold_01.cc  ---------------------------
+
+
+// Test spherical manifold on hyper shells.
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+
+// all include files you need here
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+// Helper function
+template <int dim, int spacedim>
+void test(unsigned int ref=1)
+{
+  deallog << "Testing dim " << dim 
+         << ", spacedim " << spacedim << std::endl;
+
+  SphericalManifold<dim,spacedim> manifold;
+  
+  Triangulation<dim,spacedim> tria;
+  GridGenerator::hyper_shell (tria, Point<spacedim>(), .3, .6);
+
+  for(auto cell = tria.begin_active(); cell != tria.end(); ++cell) {
+    cell->set_all_manifold_ids(1);
+  }
+  
+  tria.set_manifold(1, manifold);
+  tria.refine_global(2);
+  
+  GridOut gridout;
+  gridout.write_msh(tria, deallog.get_file_stream());
+
+  char fname[50];
+  fprints(fname, "mehs_%d_%d.msh", dim, spacedim);
+  std::ofstream of(fname);
+  gridout.write_msh(tria, of);
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  
+  test<2,2>();
+  test<3,3>();
+  
+  return 0;
+}
+
diff --git a/tests/manifold/spherical_manifold_01.output b/tests/manifold/spherical_manifold_01.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+
diff --git a/tests/manifold/spherical_manifold_02.cc b/tests/manifold/spherical_manifold_02.cc
new file mode 100644 (file)
index 0000000..5b14206
--- /dev/null
@@ -0,0 +1,74 @@
+//----------------------------  spherical_manifold_01.cc  ---------------------------
+//    Copyright (C) 2011, 2013 by the mathLab team.
+//
+//    This file is subject to LGPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  spherical_manifold_02.cc  ---------------------------
+
+
+// Test that the flat manifold does what it should on a sphere. 
+
+#include "../tests.h"
+
+#include <fstream>
+#include <base/logstream.h>
+
+
+// all include files you need here
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+// Helper function
+template <int dim, int spacedim>
+void test(unsigned int ref=1)
+{
+  SphericalManifold<dim,spacedim> manifold;
+  
+  Triangulation<dim,spacedim> tria;
+  GridGenerator::hyper_ball (tria);
+  
+  typename Triangulation<dim,spacedim>::active_cell_iterator cell;
+  
+  for(cell = tria.begin_active(); cell != tria.end(); ++cell) 
+    cell->set_all_manifold_ids(1);
+  
+  for(cell = tria.begin_active(); cell != tria.end(); ++cell) {
+    if(cell->center().distance(Point<spacedim>()) < 1e-10)
+      cell->set_all_manifold_ids(0);
+  }
+  
+  tria.set_manifold(1, manifold);
+  tria.refine_global(2);
+  
+  GridOut gridout;
+  gridout.write_msh(tria, deallog.get_file_stream());
+
+
+  char fname[50];
+  fprints(fname, "mehs_%d_%d.msh", dim, spacedim);
+  std::ofstream of(fname);
+  gridout.write_msh(tria, of);
+
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  
+  test<2,2>();
+  test<3,3>();
+  
+  return 0;
+}
+
diff --git a/tests/manifold/spherical_manifold_02.output b/tests/manifold/spherical_manifold_02.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+
diff --git a/tests/manifold/spherical_manifold_03.cc b/tests/manifold/spherical_manifold_03.cc
new file mode 100644 (file)
index 0000000..77eeb30
--- /dev/null
@@ -0,0 +1,81 @@
+//----------------------------  spherical_manifold_03.cc  ---------------------------
+//    Copyright (C) 2011, 2013 by the mathLab team.
+//
+//    This file is subject to LGPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  spherical_manifold_03.cc  ---------------------------
+
+
+// Test the push_forward and pull_back mechanisms
+
+#include "../tests.h"
+#include <fstream>
+#include <base/logstream.h>
+
+
+// all include files you need here
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
+#include <deal.II/grid/manifold_lib.h>
+#include <deal.II/grid/grid_out.h>
+
+// Helper function
+template <int dim, int spacedim>
+void test(unsigned int ref=1)
+{
+  deallog << "Testing dim " << dim
+         << ", spacedim " << spacedim << std::endl;
+
+  SphericalManifold<dim,spacedim> manifold;
+  
+  Triangulation<dim,spacedim> tria;
+  Point<spacedim> p0;
+  Point<spacedim> p1;
+  p0[0] = .5;
+  p1[0] = 1;
+  
+  if(spacedim == 2) {
+    p1[1] = 2*numbers::PI-1e-10; // theta
+  } else if(spacedim == 3) {
+    p1[1] = numbers::PI-1e-10;
+    
+    p1[2] = 2*numbers::PI-1e-10;
+  }
+
+  GridGenerator::hyper_rectangle (tria, p0, p1);
+  tria.refine_global(3);
+
+  const std::vector<Point<spacedim> > & vertices = tria.get_vertices();
+  
+  for(unsigned int i=0; i<vertices.size(); ++i) {
+    Point<spacedim> p0 = manifold.push_forward(vertices[i]);
+    Point<spacedim> p1 = manifold.pull_back(p0);
+    
+    if(p1.distance(vertices[i]) > 1e-10)
+      deallog << "ERROR! d: " << p1.distance(vertices[i]) 
+             << " - " << p1 << " != " << vertices[i] << std::endl;
+  }
+    
+  
+  
+}
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+  
+  test<2,2>();
+  test<3,3>();
+
+  return 0;
+}
+
diff --git a/tests/manifold/spherical_manifold_03.output b/tests/manifold/spherical_manifold_03.output
new file mode 100644 (file)
index 0000000..8b13789
--- /dev/null
@@ -0,0 +1 @@
+

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.