From ba1d6b53458f9e7ff595fc0587132c16ab03e846 Mon Sep 17 00:00:00 2001 From: bangerth Date: Tue, 15 Jul 2014 21:53:17 +0000 Subject: [PATCH] Update it a bit. git-svn-id: https://svn.dealii.org/trunk@33181 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-53/step-53.cc | 253 ++++++++++++++++++++++++++-- 1 file changed, 243 insertions(+), 10 deletions(-) diff --git a/deal.II/examples/step-53/step-53.cc b/deal.II/examples/step-53/step-53.cc index c4602d6720..ed19329fa0 100644 --- a/deal.II/examples/step-53/step-53.cc +++ b/deal.II/examples/step-53/step-53.cc @@ -24,12 +24,241 @@ #include #include +#include + #include using namespace dealii; +template +class SphereGeometry : public Boundary +{ +public: + SphereGeometry (const Point ¢er); + virtual + Point + get_new_point_on_line (const typename Triangulation::line_iterator &line) const; + + virtual + Point + get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const; +private: + + std_cxx1x::array pull_back (const Point &p) const; + Point push_forward (const std_cxx1x::array &preimage) const; + + template + static std_cxx1x::array average (const std_cxx1x::array (&array)[N]); + + const Point center; +}; + + +template +SphereGeometry::SphereGeometry (const Point ¢er) +: +center (center) +{} + + +template <> +std_cxx1x::array +SphereGeometry<2>::pull_back (const Point<2> &p) const +{ + const Point<2> relative_point = p - center; + + const double r = relative_point.norm(); + const double phi = std::atan2 (relative_point[1], relative_point[0]); + + std_cxx1x::array result; + result[0] = r; + result[1] = phi; + + return result; +} + + +template <> +Point<2> +SphereGeometry<2>::push_forward (const std_cxx1x::array &preimage) const +{ + const Point<2> relative_point = preimage[0] * Point<2>(std::cos(preimage[1]), std::sin(preimage[1])); + + return relative_point + center; +} + + + +template <> +std_cxx1x::array +SphereGeometry<3>::pull_back (const Point<3> &p) const +{ + const Point<3> relative_point = p - center; + + const double r = relative_point.norm(); + const double phi = std::atan2 (relative_point[1], relative_point[0]); + const double theta = std::atan2 (relative_point[2], std::sqrt (relative_point[0]*relative_point[0] + + relative_point[1]*relative_point[1])); + + std_cxx1x::array result; + result[0] = r; + result[1] = phi; + result[2] = theta; + + return result; +} + + +template <> +Point<3> +SphereGeometry<3>::push_forward (const std_cxx1x::array &preimage) const +{ + const Point<3> relative_point = (preimage[0] * + Point<3>(std::cos(preimage[1]), + std::sin(preimage[1]), + std::cos(preimage[2]))); + + return relative_point + center; +} + + + +template <> +template +std_cxx1x::array +SphereGeometry<2>::average (const std_cxx1x::array (&array)[N]) +{ + std_cxx1x::array result; + + // average the radii first. this is uncritical + { + result[0] = 0; + for (unsigned int i=0; i 1e-10) + { + const double angle = array[i][1]; + if (angle - array[0][1] > numbers::PI) + result[1] += angle-2*numbers::PI; + else if (angle - array[0][1] < -numbers::PI) + result[1] += angle+2*numbers::PI; + else + result[1] += angle; + } + else + origin_is_one_point = true; + + if (origin_is_one_point == false) + result[1] /= N; + else + result[1] /= (N-1); + } + + return result; +} + + + +template <> +template +std_cxx1x::array +SphereGeometry<3>::average (const std_cxx1x::array (&array)[N]) +{ + std_cxx1x::array result; + + // average the radii first. this is uncritical + { + result[0] = 0; + for (unsigned int i=0; i 1e-10) + { + const double angle = array[i][1]; + if (angle - array[0][1] > numbers::PI) + result[1] += angle-2*numbers::PI; + else if (angle - array[0][1] < -numbers::PI) + result[1] += angle+2*numbers::PI; + else + result[1] += angle; + } + else + origin_is_one_point = true; + + if (origin_is_one_point == false) + result[1] /= N; + else + result[1] /= (N-1); + } + + // finally for the polar angle. the difficulty here is that we have, for + // example, two points at (r,phi,theta) and (r,phi+pi,\pm pi/2), then we want + // to average these to (r,*,pi) where the equatorial angle does not matter + { + //??? not sure how exactly to do this + } + + + return result; +} + + + + +template +Point +SphereGeometry:: +get_new_point_on_line (const typename Triangulation::line_iterator &line) const +{ + std_cxx1x::array preimages[2] = { pull_back (line->vertex(0)), + pull_back (line->vertex(1)) }; + + return push_forward(average (preimages)); +} + + +template +Point +SphereGeometry:: +get_new_point_on_quad (const typename Triangulation::quad_iterator &quad) const +{ + std_cxx1x::array preimages[4] = { pull_back (quad->vertex(0)), + pull_back (quad->vertex(1)), + pull_back (quad->vertex(2)), + pull_back (quad->vertex(3)) }; + + return push_forward(average(preimages)); +} + + + template void make_grid () @@ -37,23 +266,27 @@ void make_grid () Point center; for (unsigned int i=0; i boundary(center, .25*std::sqrt((double)dim)); - Triangulation triangulation; + SphereGeometry geometry(center); + Triangulation triangulation; GridGenerator::hyper_cube (triangulation); triangulation.refine_global(1); for (typename Triangulation::active_cell_iterator cell=triangulation.begin_active(); cell!=triangulation.end(); ++cell) - for (unsigned int f=0; f::faces_per_cell; ++f) - if (cell->face(f)->center().distance(center)< radius) - cell->face(f)->set_all_manifold_ids(1); + { + if (cell->center().distance(center)< radius) + cell->set_manifold_id(1); - triangulation.set_manifold(1,boundary); - triangulation.refine_global(3); + for (unsigned int f=0; f::faces_per_cell; ++f) + if (cell->face(f)->center().distance(center)< radius) + cell->face(f)->set_all_manifold_ids(1); + } + + triangulation.set_manifold(1,geometry); + triangulation.refine_global(1); const std::string filename = "mesh-" + Utilities::int_to_string(dim) + "d.vtk"; std::ofstream out (filename.c_str()); @@ -71,5 +304,5 @@ void make_grid () int main () { make_grid<2> (); - make_grid<3> (); +// make_grid<3> (); } -- 2.39.5