Refining a hyper_ball with a spherical_manifold produces nan values

[Alexander Gary Zimmerman]

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?

[Martin Kronbichler]

Dear Alex,

you cannot set a SphericalManifold for the innermost part of a
hyper_ball because the transformation degenerates for radius=0. You need
to restrict the manifold to only the boundary or only the outer part of
the hyper_ball.

Best,
Martin

> 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 mailing list is the correct place for this kind of questions.

[luca.heltai]

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.

> --
> 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 dealii+un...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Alexander Gary Zimmerman]

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 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.

[Wolfgang Bangerth]

On 04/26/2017 02:35 AM, Alexander Gary Zimmerman wrote:
>
> 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.

Indeed. Can you open a github issue at
https://github.com/dealii/dealii/issues/new about this? Of course, we
would also gladly accept a patch to the documentation!

> 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.

That would mean that you have a hole in your domain -- that's probably
not what you want. But I think in reality, you want to read the
discussion here:

https://www.dealii.org/8.5.0/doxygen/deal.II/step_6.html#Abettermesh

Best
W.



--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/