Dear all,

attached is a minimal example which demonstrates that the function 
"extract_surface_mesh()" in the GridGenerator namespace 
does not extract refined cells. 
That is, it seems to work only for globally refined meshes. 
There is a comment in the implementation of this function:
// Create boundary mesh and mapping
// from only level(0) cells of the volume mesh.

Since the documentation states nothing about the usage of this function for 
adaptively
refined meshes (like it is done in flatten_triangulation, for instance), 
I am not sure whether this is desired behavior or not.

In my case, I probably do not need the full power of this function. 
All I want to do is to "copy" the solution values on a surface of the 
volume mesh 
to a surface mesh, and to associate a DoFHandler object to the surface mesh
that allows me to do further postprocessing.
So I do not care about the boundary and manifold id as well
as the return type, however, my volume mesh may contain hanging nodes.



Best regards,
Simon

-- 
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+unsubscr...@googlegroups.com.
To view this discussion on the web visit 
https://groups.google.com/d/msgid/dealii/f3252992-4761-41b1-a691-05c5265f89bfn%40googlegroups.com.
#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/grid_tools.h>

#include <deal.II/base/geometry_info.h>

#include <set>

using namespace dealii;

int main()
{

  Triangulation<3> tria, triaTmp;
  
  dealii::Point<3> center;
  double radius = 2.5;
  unsigned int boundaryID = 2;

  GridGenerator::hyper_cube_with_cylindrical_hole(triaTmp,
                                                    radius,
                                                    7.5,
                                                    2.0,
                                                    1,
                                                    false);
  triaTmp.refine_global(2);

  // symmetry in global Z
  std::set<typename dealii::Triangulation<3>::active_cell_iterator> cellsToRemove;
  for (const auto &cell : triaTmp.active_cell_iterators())
    if ( cell->center()[2] > 1.0)
      cellsToRemove.insert(cell);
  dealii::GridGenerator::create_triangulation_with_removed_cells(triaTmp,
                                                                  cellsToRemove,
                                                                  tria);
  // one local refinement along hole
  for (const auto &cell : tria.active_cell_iterators())
    for (unsigned int f = 0; f < dealii::GeometryInfo<3>::faces_per_cell; ++f)
      {
        dealii::Point<3> faceCenter(cell->face(f)->center());
        faceCenter[2] = 0.0;
        if ((std::abs(center.distance(faceCenter)) < radius))
          cell->set_refine_flag();
      }
  tria.execute_coarsening_and_refinement();

  // set boundary ID for surface mesh
  for (const auto &cell : tria.active_cell_iterators())
    for (unsigned int f = 0; f < dealii::GeometryInfo<3>::faces_per_cell; ++f)
      if(std::abs(cell->face(f)->center()[2]) < 1e-5 && cell->face(f)->at_boundary())
        cell->face(f)->set_boundary_id(boundaryID);        

  // extract surface mesh
  dealii::Triangulation<2, 3> tria2d;
  dealii::GridGenerator::extract_boundary_mesh<dealii::Triangulation, 3, 3>(
    tria,
    tria2d,
    std::set<dealii::types::boundary_id>{boundaryID});

  // output the volume mesh
  std::ofstream             file("tria.vtk");
  dealii::GridOut           gridOut;
  dealii::GridOutFlags::Vtk vtkFlags(true, true, true, false);
  gridOut.set_flags(vtkFlags);
  gridOut.write_vtk(tria, file);

  // and the surface mesh
  std::ofstream             file2d("tria2d.vtk");
  gridOut.write_vtk(tria2d, file2d);

  return 0;
}

Reply via email to