--- /dev/null
+
+/*
+ Code for testing the function
+ GridTools::extract_boundary_mesh (...).
+ We test that the order of cells and the orientation
+ of the vertices is consistent between the two meshes.
+*/
+
+
+#include "../tests.h"
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_out.h>
+#include <grid/grid_tools.h>
+
+using namespace std;
+
+
+template <int s_dim, int spacedim>
+bool test_vertices_orientation(const Triangulation<s_dim,spacedim> &boundary_mesh,
+ map< typename Triangulation<s_dim,spacedim>::cell_iterator,
+ typename Triangulation<s_dim+1,spacedim>::face_iterator >
+ &surface_to_volume_mapping,
+ const int verbosity = 1)
+{
+ typename Triangulation<s_dim,spacedim>::active_cell_iterator
+ cell = boundary_mesh.begin_active(),
+ endc = boundary_mesh.end();
+ typename Triangulation<s_dim+1,spacedim>::face_iterator face;
+
+ bool success = true;
+
+ if (verbosity>1)
+ deallog << "Vol faces" << "\t" << "Surf cell" <<
+ "\t" << "Distance" <<endl;
+
+ for (; cell!=endc; ++cell){
+
+ face = surface_to_volume_mapping [cell];
+
+ for (unsigned int k=0; k<GeometryInfo<s_dim>::vertices_per_cell; ++k)
+ {
+ Point<spacedim> diff(face->vertex(k));
+ diff -= cell->vertex(k);
+ if (verbosity>1){
+ deallog << face->vertex(k) << "\t\t";
+ deallog << cell->vertex(k) << "\t\t\t" << diff.square() << endl;
+ }
+ if (diff.square()>0){
+ success = false;
+ break;
+ }
+ }
+ if (verbosity>1) deallog << endl;
+ }
+ return success;
+}
+
+template <int dim, int spacedim>
+void save_mesh(const Triangulation<dim,spacedim>& tria)
+{
+ GridOut grid_out;
+ grid_out.write_ucd (tria, deallog.get_file_stream());
+}
+
+
+int main ()
+{
+
+ ofstream logfile("extract_boundary_mesh_00/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+ { // Extract the whole boundary of a hyper-cube
+ const int dim = 2;
+
+ deallog << "Testing hyper_cube in dim: " << dim << "..."<< endl;
+ map< Triangulation<dim-1,dim>::cell_iterator,
+ Triangulation<dim,dim>::face_iterator>
+ surface_to_volume_mapping;
+
+ Triangulation<dim> volume_mesh;
+ GridGenerator::hyper_cube(volume_mesh);
+ volume_mesh.begin_active()->face(0)->set_boundary_indicator(1);
+ volume_mesh.refine_global (1);
+
+ save_mesh(volume_mesh);
+
+ Triangulation<dim-1,dim> boundary_mesh;
+
+
+ GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh,
+ surface_to_volume_mapping);
+
+ if (test_vertices_orientation(boundary_mesh, surface_to_volume_mapping))
+ deallog << "Passed.";
+ else
+ deallog << "Failed.";
+ deallog << endl;
+
+ save_mesh(boundary_mesh);
+
+ }
+
+ { // Extract the whole boundary of a hyper-cube
+ const int dim = 3;
+
+ deallog << "Testing hyper_cube in dim: " << dim << "..."<< endl;
+ map< Triangulation<dim-1,dim>::cell_iterator,
+ Triangulation<dim,dim>::face_iterator>
+ surface_to_volume_mapping;
+
+ Triangulation<dim> volume_mesh;
+ GridGenerator::hyper_cube(volume_mesh);
+ volume_mesh.begin_active()->face(0)->set_boundary_indicator(1);
+ volume_mesh.refine_global (1);
+
+ save_mesh(volume_mesh);
+
+ Triangulation<dim-1,dim> boundary_mesh;
+
+
+ GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh,
+ surface_to_volume_mapping);
+
+ if (test_vertices_orientation(boundary_mesh, surface_to_volume_mapping))
+ deallog << "Passed.";
+ else
+ deallog << "Failed.";
+ deallog << endl;
+
+ save_mesh(boundary_mesh);
+
+ }
+
+
+ { // Extract a piece of the boundary of a hyper-cube
+
+ const int dim = 3;
+ deallog << "Testing hyper_cube in dim: " << dim << "..."<< endl;
+
+ map< Triangulation<dim-1,dim>::cell_iterator,
+ Triangulation<dim,dim>::face_iterator>
+ surface_to_volume_mapping;
+
+ Triangulation<dim> volume_mesh;
+ GridGenerator::hyper_cube(volume_mesh);
+ volume_mesh.begin_active()->face(0)->set_boundary_indicator(1);
+ volume_mesh.refine_global (1);
+
+ save_mesh(volume_mesh);
+
+ Triangulation<dim-1,dim> boundary_mesh;
+ set<unsigned char> boundary_ids;
+ boundary_ids.insert(0);
+
+ GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh,
+ surface_to_volume_mapping,
+ boundary_ids);
+
+ if (test_vertices_orientation(boundary_mesh, surface_to_volume_mapping))
+ deallog << "Passed.";
+ else
+ deallog << "Failed.";
+ deallog << endl;
+
+ save_mesh(boundary_mesh);
+
+ }
+
+
+
+
+
+ return 0;
+ }
--- /dev/null
+
+DEAL::Testing hyper_cube in dim: 2...
+9 4 0 0 0
+1 0.00000 0.00000 0
+2 1.00000 0.00000 0
+3 0.00000 1.00000 0
+4 1.00000 1.00000 0
+5 0.500000 0.00000 0
+6 0.00000 0.500000 0
+7 1.00000 0.500000 0
+8 0.500000 1.00000 0
+9 0.500000 0.500000 0
+1 0 quad 1 5 9 6
+2 0 quad 5 2 7 9
+3 0 quad 6 9 8 3
+4 0 quad 9 7 4 8
+DEAL::Passed.
+8 8 0 0 0
+1 0.00000 0.00000 0
+2 0.00000 1.00000 0
+3 1.00000 0.00000 0
+4 1.00000 1.00000 0
+5 0.00000 0.500000 0
+6 1.00000 0.500000 0
+7 0.500000 0.00000 0
+8 0.500000 1.00000 0
+1 0 line 1 5
+2 0 line 5 2
+3 0 line 3 6
+4 0 line 6 4
+5 0 line 1 7
+6 0 line 7 3
+7 0 line 2 8
+8 0 line 8 4
+DEAL::Testing hyper_cube in dim: 3...
+27 8 0 0 0
+1 0.00000 0.00000 0.00000
+2 1.00000 0.00000 0.00000
+3 0.00000 1.00000 0.00000
+4 1.00000 1.00000 0.00000
+5 0.00000 0.00000 1.00000
+6 1.00000 0.00000 1.00000
+7 0.00000 1.00000 1.00000
+8 1.00000 1.00000 1.00000
+9 0.500000 0.00000 0.00000
+10 0.00000 0.500000 0.00000
+11 0.00000 0.00000 0.500000
+12 1.00000 0.500000 0.00000
+13 1.00000 0.00000 0.500000
+14 0.500000 1.00000 0.00000
+15 0.00000 1.00000 0.500000
+16 1.00000 1.00000 0.500000
+17 0.500000 0.00000 1.00000
+18 0.00000 0.500000 1.00000
+19 1.00000 0.500000 1.00000
+20 0.500000 1.00000 1.00000
+21 0.500000 0.00000 0.500000
+22 0.500000 0.500000 0.00000
+23 0.00000 0.500000 0.500000
+24 1.00000 0.500000 0.500000
+25 0.500000 1.00000 0.500000
+26 0.500000 0.500000 1.00000
+27 0.500000 0.500000 0.500000
+1 0 hex 1 9 21 11 10 22 27 23
+2 0 hex 9 2 13 21 22 12 24 27
+3 0 hex 10 22 27 23 3 14 25 15
+4 0 hex 22 12 24 27 14 4 16 25
+5 0 hex 11 21 17 5 23 27 26 18
+6 0 hex 21 13 6 17 27 24 19 26
+7 0 hex 23 27 26 18 15 25 20 7
+8 0 hex 27 24 19 26 25 16 8 20
+DEAL::Passed.
+26 24 0 0 0
+1 0.00000 0.00000 0.00000
+2 0.00000 1.00000 0.00000
+3 0.00000 0.00000 1.00000
+4 0.00000 1.00000 1.00000
+5 1.00000 0.00000 0.00000
+6 1.00000 1.00000 0.00000
+7 1.00000 0.00000 1.00000
+8 1.00000 1.00000 1.00000
+9 0.00000 0.500000 0.00000
+10 0.00000 0.00000 0.500000
+11 0.500000 0.00000 0.00000
+12 0.00000 1.00000 0.500000
+13 0.500000 1.00000 0.00000
+14 0.00000 0.500000 1.00000
+15 0.500000 0.00000 1.00000
+16 0.500000 1.00000 1.00000
+17 1.00000 0.500000 0.00000
+18 1.00000 0.00000 0.500000
+19 1.00000 1.00000 0.500000
+20 1.00000 0.500000 1.00000
+21 0.00000 0.500000 0.500000
+22 1.00000 0.500000 0.500000
+23 0.500000 0.00000 0.500000
+24 0.500000 1.00000 0.500000
+25 0.500000 0.500000 0.00000
+26 0.500000 0.500000 1.00000
+1 0 quad 1 9 21 10
+2 0 quad 9 2 12 21
+3 0 quad 10 21 14 3
+4 0 quad 21 12 4 14
+5 0 quad 5 17 22 18
+6 0 quad 17 6 19 22
+7 0 quad 18 22 20 7
+8 0 quad 22 19 8 20
+9 0 quad 1 10 23 11
+10 0 quad 10 3 15 23
+11 0 quad 11 23 18 5
+12 0 quad 23 15 7 18
+13 0 quad 2 12 24 13
+14 0 quad 12 4 16 24
+15 0 quad 13 24 19 6
+16 0 quad 24 16 8 19
+17 0 quad 1 11 25 9
+18 0 quad 11 5 17 25
+19 0 quad 9 25 13 2
+20 0 quad 25 17 6 13
+21 0 quad 3 15 26 14
+22 0 quad 15 7 20 26
+23 0 quad 14 26 16 4
+24 0 quad 26 20 8 16
+DEAL::Testing hyper_cube in dim: 3...
+27 8 0 0 0
+1 0.00000 0.00000 0.00000
+2 1.00000 0.00000 0.00000
+3 0.00000 1.00000 0.00000
+4 1.00000 1.00000 0.00000
+5 0.00000 0.00000 1.00000
+6 1.00000 0.00000 1.00000
+7 0.00000 1.00000 1.00000
+8 1.00000 1.00000 1.00000
+9 0.500000 0.00000 0.00000
+10 0.00000 0.500000 0.00000
+11 0.00000 0.00000 0.500000
+12 1.00000 0.500000 0.00000
+13 1.00000 0.00000 0.500000
+14 0.500000 1.00000 0.00000
+15 0.00000 1.00000 0.500000
+16 1.00000 1.00000 0.500000
+17 0.500000 0.00000 1.00000
+18 0.00000 0.500000 1.00000
+19 1.00000 0.500000 1.00000
+20 0.500000 1.00000 1.00000
+21 0.500000 0.00000 0.500000
+22 0.500000 0.500000 0.00000
+23 0.00000 0.500000 0.500000
+24 1.00000 0.500000 0.500000
+25 0.500000 1.00000 0.500000
+26 0.500000 0.500000 1.00000
+27 0.500000 0.500000 0.500000
+1 0 hex 1 9 21 11 10 22 27 23
+2 0 hex 9 2 13 21 22 12 24 27
+3 0 hex 10 22 27 23 3 14 25 15
+4 0 hex 22 12 24 27 14 4 16 25
+5 0 hex 11 21 17 5 23 27 26 18
+6 0 hex 21 13 6 17 27 24 19 26
+7 0 hex 23 27 26 18 15 25 20 7
+8 0 hex 27 24 19 26 25 16 8 20
+DEAL::Passed.
+25 20 0 0 0
+1 1.00000 0.00000 0.00000
+2 1.00000 1.00000 0.00000
+3 1.00000 0.00000 1.00000
+4 1.00000 1.00000 1.00000
+5 0.00000 0.00000 0.00000
+6 0.00000 0.00000 1.00000
+7 0.00000 1.00000 0.00000
+8 0.00000 1.00000 1.00000
+9 1.00000 0.500000 0.00000
+10 1.00000 0.00000 0.500000
+11 1.00000 1.00000 0.500000
+12 1.00000 0.500000 1.00000
+13 0.500000 0.00000 0.00000
+14 0.00000 0.00000 0.500000
+15 0.00000 0.500000 0.00000
+16 0.500000 0.00000 1.00000
+17 0.00000 0.500000 1.00000
+18 0.500000 1.00000 0.00000
+19 0.00000 1.00000 0.500000
+20 0.500000 1.00000 1.00000
+21 1.00000 0.500000 0.500000
+22 0.500000 0.00000 0.500000
+23 0.500000 1.00000 0.500000
+24 0.500000 0.500000 0.00000
+25 0.500000 0.500000 1.00000
+1 0 quad 1 9 21 10
+2 0 quad 9 2 11 21
+3 0 quad 10 21 12 3
+4 0 quad 21 11 4 12
+5 0 quad 5 14 22 13
+6 0 quad 14 6 16 22
+7 0 quad 13 22 10 1
+8 0 quad 22 16 3 10
+9 0 quad 7 19 23 18
+10 0 quad 19 8 20 23
+11 0 quad 18 23 11 2
+12 0 quad 23 20 4 11
+13 0 quad 5 13 24 15
+14 0 quad 13 1 9 24
+15 0 quad 15 24 18 7
+16 0 quad 24 9 2 18
+17 0 quad 6 16 25 17
+18 0 quad 16 3 12 25
+19 0 quad 17 25 20 8
+20 0 quad 25 12 4 20
--- /dev/null
+
+/*
+ Code for testing the function
+ GridTools::extract_boundary_mesh (...).
+ We test that the order of cells and the orientation
+ of the vertices is consistent between the two meshes.
+
+ It's faling when boundary object is attached and refinement
+ for a ball in 3D
+*/
+
+#include "../tests.h"
+
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/grid_out.h>
+#include <grid/grid_tools.h>
+
+using namespace std;
+
+
+template <int s_dim, int spacedim>
+bool test_vertices_orientation(const Triangulation<s_dim,spacedim> &bdry_mesh,
+ map< typename Triangulation<s_dim,spacedim>::cell_iterator,
+ typename Triangulation<s_dim+1,spacedim>::face_iterator >
+ &surface_to_volume_mapping,
+ const int verbosity = 1)
+{
+ typename Triangulation<s_dim,spacedim>::active_cell_iterator
+ cell = bdry_mesh.begin_active(),
+ endc = bdry_mesh.end();
+ typename Triangulation<s_dim+1,spacedim>::face_iterator face;
+
+ bool success = true;
+
+ if (verbosity>1){
+ deallog << "The last column should be 0" <<endl;
+ deallog << "Vol faces" << "\t\t" << "Surf cell" <<
+ "\t\t" << "Distance" <<endl<<endl;
+ }
+
+ for (; cell!=endc; ++cell){
+
+ face = surface_to_volume_mapping [cell];
+
+ for (unsigned int k=0; k<GeometryInfo<s_dim>::vertices_per_cell; ++k)
+ {
+ Point<spacedim> diff(face->vertex(k));
+ diff -= cell->vertex(k);
+ if (verbosity>1){
+ deallog << face->vertex(k) << "\t\t";
+ deallog << cell->vertex(k) << "\t\t\t" << diff.square() << endl;
+ }
+ if (diff.square()>0){
+ success = false;
+ break;
+ }
+ }
+ if (verbosity>1) deallog << endl;
+ }
+ return success;
+}
+
+template <int dim, int spacedim>
+void save_mesh(const Triangulation<dim,spacedim>& tria)
+{
+ GridOut grid_out;
+ grid_out.write_ucd (tria, deallog.get_file_stream());
+}
+
+
+int main ()
+{
+
+ ofstream logfile("extract_boundary_mesh_01/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+
+
+ { // Extract the boundary of a hyper-sphere
+
+ const int dim = 3;
+ deallog << "Testing hyper_cube in dim: " << dim << "..."<< endl;
+
+ map< Triangulation<dim-1,dim>::cell_iterator,
+ Triangulation<dim,dim>::face_iterator>
+ surface_to_volume_mapping;
+ const HyperBallBoundary<dim> boundary_description;
+ Triangulation<dim> volume_mesh;
+ GridGenerator::hyper_ball(volume_mesh);
+ volume_mesh.set_boundary (0, boundary_description);
+
+ volume_mesh.refine_global (1);
+
+ // save_mesh(volume_mesh);
+
+ Triangulation<dim-1,dim> bdry_mesh;
+
+ GridTools::extract_boundary_mesh (volume_mesh, bdry_mesh, surface_to_volume_mapping);
+
+ if (test_vertices_orientation(bdry_mesh, surface_to_volume_mapping,2))
+ deallog << "Passed.";
+ else
+ deallog << "Failed.";
+ deallog << endl;
+
+ // save_mesh(bdry_mesh);
+
+ }
+
+
+ return 0;
+ }