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].
For more options, visit https://groups.google.com/d/optout.

Reply via email to