Thanks, Luca. The behavior is obvious to me now. As for the documentation: I indeed missed that detail from SphericalManifold; but* there is a very misleading statement in the hyper_ball documentation:* "You should attach a SphericalManifold <http://dealii.org/8.5.0/doxygen/deal.II/classSphericalManifold.html> to the cells and faces for correct placement of vertices upon refinement and to be able to use higher order mappings." I would guess that this could lead many people to (at least initially) make my mistake, so maybe a small edit to note the pitfall could be made here.
So if I want to have a solid sphere domain with a spherical manifold throughout it's entirety, I understand that my best option is to use a hyper_shell with a very small inside radius. On Wednesday, April 26, 2017 at 10:06:00 AM UTC+2, Luca Heltai wrote: > > This is actually not a bug. A spherical manifold is undefined on the > center of the sphere/ball. You cannot attach set_all_manifold_ids of the > spherical manifold and attach a spherical manifold to it, since the central > cell (the one where the cell->center() coincides with the ball center) > would not know what to do. > > This is documented in the SphericalManfold documentation. See the phrase > "In order for this manifold to work correctly, it cannot be attached to > cells containing the center of the coordinate system. This point is a > singular point of the coordinate transformation, and there taking averages > does not make any sense.” > > Best, > > L. > > > > On 25 Apr 2017, at 13:42, Alexander Gary Zimmerman < > [email protected] <javascript:>> wrote: > > > > When I refine a Triangulation that has been initialized with > GridGenerator::hyper_ball, and has a SphericalManifold attached (and having > set_all_manifold_ids on the Triangulation), some nan (not a number) values > are produced. > > > > On the other hand, I get the expected behavior when instead using a > GridGenerator::hyper_shell. > > > > I've documented the issue on my repo here: > https://github.com/alexanderzimmerman/nsb-pcm/issues/21 > > > > I'm using deal.II-v.8.5.0. I haven't tried reproducing this on any other > versions. > > > > These are in the repo, but here's the test that passes: > > > > #include <iostream> > > #include <fstream> > > #include <string> > > > > #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/grid_refinement.h> > > #include <deal.II/grid/grid_out.h> > > #include <deal.II/grid/manifold_lib.h> > > #include <deal.II/base/numbers.h> > > #include <deal.II/base/signaling_nan.h> > > > > using namespace dealii; > > > > const unsigned int dim = 2; > > > > > > void validate_vertices(Triangulation<dim> &tria) > > { > > Triangulation<dim>::cell_iterator cell; > > > > unsigned int nv = GeometryInfo<dim>::vertices_per_cell; > > > > for (cell = tria.begin(); cell != tria.end(); ++cell) > > { > > for (unsigned int v = 0; v < nv; ++v) > > { > > for (unsigned int i = 0; i < dim; ++i) > > { > > assert(!numbers::is_nan(cell->vertex(v)[i])); > > } > > > > } > > > > } > > > > } > > > > > > int main(int /*argc*/, char** /*argv*/) > > { > > Triangulation<dim> triangulation; > > > > GridGenerator::hyper_shell(triangulation, Point<dim>(), 0.1, 0.5); > > > > const unsigned int manifold_id = 0; > > > > triangulation.set_all_manifold_ids(manifold_id); > > > > SphericalManifold<dim> spherical_manifold; > > > > triangulation.set_manifold(manifold_id, spherical_manifold); > > > > GridOut grid_writer; > > > > > > for (unsigned int r = 0; r < 2; ++r) > > { > > std::string file_name = "grid_"+std::to_string(r)+".vtk"; > > > > std::ofstream file(file_name); > > > > grid_writer.write_vtk(triangulation, file); > > > > validate_vertices(triangulation); > > > > triangulation.refine_global(1); > > } > > > > triangulation.set_manifold(0); > > > > return 0; > > } > > > > > > and here's the test that fails: > > > > #include <iostream> > > #include <fstream> > > #include <string> > > > > #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/grid_refinement.h> > > #include <deal.II/grid/grid_out.h> > > #include <deal.II/grid/manifold_lib.h> > > #include <deal.II/base/numbers.h> > > #include <deal.II/base/signaling_nan.h> > > > > > > using namespace dealii; > > > > const unsigned int dim = 2; > > > > > > void validate_vertices(Triangulation<dim> &tria) > > { > > Triangulation<dim>::cell_iterator cell; > > > > unsigned int nv = GeometryInfo<dim>::vertices_per_cell; > > > > for (cell = tria.begin(); cell != tria.end(); ++cell) > > { > > for (unsigned int v = 0; v < nv; ++v) > > { > > for (unsigned int i = 0; i < dim; ++i) > > { > > assert(!numbers::is_nan(cell->vertex(v)[i])); > > } > > > > } > > > > } > > > > } > > > > > > int main(int /*argc*/, char** /*argv*/) > > { > > Triangulation<dim> triangulation; > > > > GridGenerator::hyper_ball(triangulation); > > > > const unsigned int manifold_id = 0; > > > > triangulation.set_all_manifold_ids(manifold_id); > > > > SphericalManifold<dim> spherical_manifold; > > > > triangulation.set_manifold(manifold_id, spherical_manifold); > > > > GridOut grid_writer; > > > > > > for (unsigned int r = 0; r < 2; ++r) > > { > > std::string file_name = "grid_"+std::to_string(r)+".vtk"; > > > > std::ofstream file(file_name); > > > > grid_writer.write_vtk(triangulation, file); > > > > validate_vertices(triangulation); > > > > triangulation.refine_global(1); > > } > > > > triangulation.set_manifold(0); > > > > return 0; > > } > > > > > > Any ideas? > > > > Thanks, > > > > Alex > > > > P.S. for future reference, if I have an issue formulated like this, > should I raise it on the deal.II repo or here in the mailing list? > > > > -- > > The deal.II project is located at http://www.dealii.org/ > > For mailing list/forum options, see > https://groups.google.com/d/forum/dealii?hl=en > > --- > > You received this message because you are subscribed to the Google > Groups "deal.II User Group" group. > > To unsubscribe from this group and stop receiving emails from it, send > an email to [email protected] <javascript:>. > > For more options, visit https://groups.google.com/d/optout. > > -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. For more options, visit https://groups.google.com/d/optout.
