Fixed compilation bug on Mac.
Modified tests for spherical manifold.
--- /dev/null
+// ---------------------------------------------------------------------
+// $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
grid_tools.cc
intergrid_map.cc
manifold.cc
+ manifold_lib.cc
persistent_tria.cc
tria_accessor.cc
tria_boundary.cc
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;
--- /dev/null
+// ---------------------------------------------------------------------
+// $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
--- /dev/null
+// ---------------------------------------------------------------------
+// $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>;
+ }
+
+
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;
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.
--- /dev/null
+//---------------------------- 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;
+}
+
--- /dev/null
+//---------------------------- 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;
+}
+
--- /dev/null
+//---------------------------- 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;
+}
+
projects / dealii.git / commitdiff