From 36c03228a2b65768e7dd25428054c2870ac5439c Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Fri, 1 Aug 2014 11:08:04 +0200 Subject: [PATCH] Added buggy spherical manifold. Fixed compilation bug on Mac. Modified tests for spherical manifold. --- include/deal.II/grid/manifold_lib.h | 87 +++++++++++++++ source/grid/CMakeLists.txt | 1 + source/grid/manifold.cc | 2 +- source/grid/manifold_lib.cc | 113 ++++++++++++++++++++ source/grid/manifold_lib.inst.in | 25 +++++ source/grid/tria_boundary_lib.cc | 4 +- tests/manifold/spherical_manifold_01.cc | 68 ++++++++++++ tests/manifold/spherical_manifold_01.output | 1 + tests/manifold/spherical_manifold_02.cc | 74 +++++++++++++ tests/manifold/spherical_manifold_02.output | 1 + tests/manifold/spherical_manifold_03.cc | 81 ++++++++++++++ tests/manifold/spherical_manifold_03.output | 1 + 12 files changed, 455 insertions(+), 3 deletions(-) create mode 100644 include/deal.II/grid/manifold_lib.h create mode 100644 source/grid/manifold_lib.cc create mode 100644 source/grid/manifold_lib.inst.in create mode 100644 tests/manifold/spherical_manifold_01.cc create mode 100644 tests/manifold/spherical_manifold_01.output create mode 100644 tests/manifold/spherical_manifold_02.cc create mode 100644 tests/manifold/spherical_manifold_02.output create mode 100644 tests/manifold/spherical_manifold_03.cc create mode 100644 tests/manifold/spherical_manifold_03.output diff --git a/include/deal.II/grid/manifold_lib.h b/include/deal.II/grid/manifold_lib.h new file mode 100644 index 0000000000..6a1f7ee63a --- /dev/null +++ b/include/deal.II/grid/manifold_lib.h @@ -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 +#include + +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, 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 +class SphericalManifold : public ManifoldChart +{ +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 center = Point()); + + /** + * Pull back the given point from the Euclidean space. Will return + * the polar coordinates associated with the point @p space_point. + */ + virtual Point + pull_back(const Point &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 + push_forward(const Point &chart_point) const; + + + /** + * The center of the spherical coordinate system. + */ + const Point center; +private: + + /** Helper function which returns the periodicity associated with + this coordinate system, according to dim, chartdim, and + spacedim. */ + static Point get_periodicity(); +}; + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/grid/CMakeLists.txt b/source/grid/CMakeLists.txt index 1b074c9523..18e30fe8c7 100644 --- a/source/grid/CMakeLists.txt +++ b/source/grid/CMakeLists.txt @@ -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 diff --git a/source/grid/manifold.cc b/source/grid/manifold.cc index 1685bbc1be..6d4dd06c81 100644 --- a/source/grid/manifold.cc +++ b/source/grid/manifold.cc @@ -64,7 +64,7 @@ namespace Manifolds { template void get_default_quadrature(const OBJECT& obj, - Quadrature &quad, bool with_laplace = false) + Quadrature &quad, bool with_laplace) { std::vector > sp; std::vector wp; diff --git a/source/grid/manifold_lib.cc b/source/grid/manifold_lib.cc new file mode 100644 index 0000000000..d410187425 --- /dev/null +++ b/source/grid/manifold_lib.cc @@ -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 +#include +#include +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +template +SphericalManifold::SphericalManifold(const Point center): + ManifoldChart(SphericalManifold::get_periodicity()), + center(center) +{ + Assert(spacedim != 1, ExcImpossibleInDim(1)); +} + + + +template +Point +SphericalManifold::get_periodicity() { + Point periodicity; + periodicity[spacedim-1] = 2*numbers::PI; // theta and phi period. + return periodicity; +} + + + +template +Point +SphericalManifold::push_forward(const Point &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 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 +Point +SphericalManifold::pull_back(const Point &space_point) const { + const Point R = space_point-center; + const double rho = R.norm(); + const double &x = R[0]; + const double &y = R[1]; + + Point 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 index 0000000000..11bf4195a9 --- /dev/null +++ b/source/grid/manifold_lib.inst.in @@ -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; + template class CylinderManifold; +// template class ConeBoundary; + } + + diff --git a/source/grid/tria_boundary_lib.cc b/source/grid/tria_boundary_lib.cc index 0043e63e3b..a2ff2f43df 100644 --- a/source/grid/tria_boundary_lib.cc +++ b/source/grid/tria_boundary_lib.cc @@ -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 axis = x_1 - x_0; // Compute the middle point of the quad. - const Point middle = StraightBoundary<3>::get_new_point_on_quad (quad); + const Point 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 index 0000000000..cf36d3b235 --- /dev/null +++ b/tests/manifold/spherical_manifold_01.cc @@ -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 +#include + + +// all include files you need here +#include +#include +#include +#include +#include +#include +#include + +// Helper function +template +void test(unsigned int ref=1) +{ + deallog << "Testing dim " << dim + << ", spacedim " << spacedim << std::endl; + + SphericalManifold manifold; + + Triangulation tria; + GridGenerator::hyper_shell (tria, Point(), .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 index 0000000000..8b13789179 --- /dev/null +++ b/tests/manifold/spherical_manifold_01.output @@ -0,0 +1 @@ + diff --git a/tests/manifold/spherical_manifold_02.cc b/tests/manifold/spherical_manifold_02.cc new file mode 100644 index 0000000000..5b14206afe --- /dev/null +++ b/tests/manifold/spherical_manifold_02.cc @@ -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 +#include + + +// all include files you need here +#include +#include +#include +#include +#include +#include +#include + +// Helper function +template +void test(unsigned int ref=1) +{ + SphericalManifold manifold; + + Triangulation tria; + GridGenerator::hyper_ball (tria); + + typename Triangulation::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()) < 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 index 0000000000..8b13789179 --- /dev/null +++ b/tests/manifold/spherical_manifold_02.output @@ -0,0 +1 @@ + diff --git a/tests/manifold/spherical_manifold_03.cc b/tests/manifold/spherical_manifold_03.cc new file mode 100644 index 0000000000..77eeb30e2c --- /dev/null +++ b/tests/manifold/spherical_manifold_03.cc @@ -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 +#include + + +// all include files you need here +#include +#include +#include +#include +#include +#include +#include + +// Helper function +template +void test(unsigned int ref=1) +{ + deallog << "Testing dim " << dim + << ", spacedim " << spacedim << std::endl; + + SphericalManifold manifold; + + Triangulation tria; + Point p0; + Point 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 > & vertices = tria.get_vertices(); + + for(unsigned int i=0; i p0 = manifold.push_forward(vertices[i]); + Point 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 index 0000000000..8b13789179 --- /dev/null +++ b/tests/manifold/spherical_manifold_03.output @@ -0,0 +1 @@ + -- 2.39.5