From c24bb1bebacc4b5479ab9d4c41be0995730afe2b Mon Sep 17 00:00:00 2001 From: heltai Date: Tue, 28 Jan 2014 10:17:57 +0000 Subject: [PATCH] Added periodicity to FlatManifold git-svn-id: https://svn.dealii.org/branches/branch_manifold_id@32331 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/include/deal.II/grid/manifold.h | 46 +++++++++++++++++--- deal.II/source/grid/manifold.cc | 29 +++++++++++-- deal.II/source/grid/tria_boundary_lib.cc | 3 +- tests/manifold/flat_manifold_01.cc | 55 ++++++++++++++++++++++++ tests/manifold/flat_manifold_01.output | 16 +++++++ 5 files changed, 137 insertions(+), 12 deletions(-) create mode 100644 tests/manifold/flat_manifold_01.cc create mode 100644 tests/manifold/flat_manifold_01.output diff --git a/deal.II/include/deal.II/grid/manifold.h b/deal.II/include/deal.II/grid/manifold.h index ed85ca3678..955962b69c 100644 --- a/deal.II/include/deal.II/grid/manifold.h +++ b/deal.II/include/deal.II/grid/manifold.h @@ -151,9 +151,22 @@ class FlatManifold: public Manifold public: /** * Default constructor. Some compilers require this for some - * reasons. + * reasons. The The optional argument can be used to specify the + * periodicity of the spacedim-dimensional manifold (one period per + * direction). A peridicity value of zero means that along that + * direction there is no peridicity. By default no periodicity is + * assumed. + * + * Periodicity affects the way a middle point is computed. It is + * assumed that if two points are more than half period distant, + * then the distance should be computed by crossing the periodicity + * boundary, i.e., the average is computed by adding a full period + * to the sum of the two. For example, if along direction 0 we have + * 2*pi periodicity, then the average of (2*pi-eps) and (eps) is not + * pi, but 2*pi (or zero), since, on a periodic manifold, these two + * points are at distance 2*eps and not (2*pi-eps) */ - FlatManifold (); + FlatManifold (const Point periodicity=Point()); /** * Let the new point be the average sum of surrounding vertices. @@ -171,7 +184,7 @@ public: virtual Point get_new_point(const std::vector > &surrounding_points, const std::vector &weights) const; - + /** * Project to FlatManifold: do nothing. Note however that this * function can be overloaded by derived classes, which will then @@ -181,6 +194,8 @@ public: */ virtual Point project_to_manifold (const Point &candidate) const; +private: + const Point periodicity; }; @@ -227,7 +242,24 @@ template class ManifoldChart: public Manifold { public: - + /** + * Constructor. The optional argument can be used to specify the + * periodicity of the dim-dimensional manifold (one period per + * direction). A peridicity value of zero means that along that + * direction there is no peridicity. By default no periodicity is + * assumed. + * + * Periodicity affects the way a middle point is computed. It is + * assumed that if two points are more than half period distant, + * then the distance should be computed by crossing the periodicity + * boundary, i.e., then the average is computed by adding a full + * period to the sum of the two. For example, if along direction 0 + * we have 2*pi periodicity, then the average of (2*pi-eps) and + * (eps) is not pi, but 2*pi (or zero), since, on the manifold, + * these two points are at distance 2*eps and not (2*pi-eps) + */ + ManifoldChart(const Point periodicity=Point()); + /** * Destructor. Does nothing here, but needs to be declared to make * it virtual. @@ -264,10 +296,12 @@ public: virtual Point push_forward(const Point &chart_point) const = 0; - private: + /** + * The sub_manifold object is used to compute the average of the + * points in the chart coordinates system. + */ FlatManifold sub_manifold; - mutable std::vector > chart_points; }; diff --git a/deal.II/source/grid/manifold.cc b/deal.II/source/grid/manifold.cc index 1d598217a8..60f6fb8f25 100644 --- a/deal.II/source/grid/manifold.cc +++ b/deal.II/source/grid/manifold.cc @@ -70,7 +70,8 @@ Manifold::normal_vector(const Point &point) const template -FlatManifold::FlatManifold () +FlatManifold::FlatManifold (const Point periodicity) : + periodicity(periodicity) {} template @@ -90,8 +91,23 @@ get_new_point (const std::vector > &surrounding_points, Point p; - for(unsigned int i=0; i dp; + bool check_period = (periodicity.norm() != 0); + for(unsigned int i=0; i 0) + dp[d] = std::abs(surrounding_points[i][d]- + surrounding_points[0][d]) > periodicity[d]/2.0 ? + periodicity[d] : 0.0; + } + p += (surrounding_points[i]+dp)*weights[i]; + } + if(check_period) + for(unsigned int d=0; d 0) + p[d] = std::fmod(p[d], periodicity[d]); + return project_to_manifold(p); } @@ -109,6 +125,11 @@ template ManifoldChart::~ManifoldChart () {} +template +ManifoldChart::ManifoldChart (const Point periodicity): + sub_manifold(periodicity) +{} + template Point @@ -117,7 +138,7 @@ get_new_point (const std::vector > &surrounding_points, const std::vector &weights) const { AssertDimension(surrounding_points.size(), weights.size()); - chart_points.resize(surrounding_points.size()); + std::vector > chart_points(surrounding_points.size()); for(unsigned int i=0; i TorusBoundary<2,3>::TorusBoundary (const double R__, const double r__) : + ManifoldChart<2,3>(Point<2>(2*numbers::PI, 2*numbers::PI)), R(R__), r(r__) { @@ -448,8 +449,6 @@ TorusBoundary::get_correct_angle(const double angle, return 2.0*numbers::PI-angle; } - - template <> Point<3> TorusBoundary<2,3>::push_forward (const Point<2> &surfP) const diff --git a/tests/manifold/flat_manifold_01.cc b/tests/manifold/flat_manifold_01.cc new file mode 100644 index 0000000000..9d463df4e4 --- /dev/null +++ b/tests/manifold/flat_manifold_01.cc @@ -0,0 +1,55 @@ +//---------------------------- manifold_id_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. +// +//---------------------------- flat_manifold_01.cc --------------------------- + + +// Test periodicity conditions using FlatManifold + +#include "../tests.h" +#include +#include + + +// all include files you need here +#include +#include + +template +void test() { + Point period; + for(unsigned int i=0; i manifold(period); + deallog << "Testing spacedim = " << spacedim << std::endl; + deallog << "Period: " << period << std::endl; + std::vector > p(2); + std::vector w(2, 0.5); + double eps = .1; + p[0] = period*eps*2; + p[1] = period*(1.0-eps); + Point middle = manifold.get_new_point(p, w); + deallog << "P(0): " << p[0] << std::endl; + deallog << "P(1): " << p[1] << std::endl; + deallog << "Middle: " << middle << std::endl; +} + +int main () +{ + std::ofstream logfile("output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-5); + + test<1>(); + test<2>(); + test<3>(); + + return 0; +} + diff --git a/tests/manifold/flat_manifold_01.output b/tests/manifold/flat_manifold_01.output new file mode 100644 index 0000000000..f6f885fc5b --- /dev/null +++ b/tests/manifold/flat_manifold_01.output @@ -0,0 +1,16 @@ + +DEAL::Testing spacedim = 1 +DEAL::Period: 1.00000 +DEAL::P(0): 0.200000 +DEAL::P(1): 0.900000 +DEAL::Middle: 0.0500000 +DEAL::Testing spacedim = 2 +DEAL::Period: 1.00000 1.00000 +DEAL::P(0): 0.200000 0.200000 +DEAL::P(1): 0.900000 0.900000 +DEAL::Middle: 0.0500000 0.0500000 +DEAL::Testing spacedim = 3 +DEAL::Period: 1.00000 1.00000 1.00000 +DEAL::P(0): 0.200000 0.200000 0.200000 +DEAL::P(1): 0.900000 0.900000 0.900000 +DEAL::Middle: 0.0500000 0.0500000 0.0500000 -- 2.39.5