Extracting a piece of the boundary of the sphere does not seem

to project some nodes properly

git-svn-id: https://svn.dealii.org/trunk@22809 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/tests/codim_one/extract_boundary_mesh_04.cc b/tests/codim_one/extract_boundary_mesh_04.cc
new file mode 100644 (file)
index 0000000..6a03470
--- /dev/null
+++ b/tests/codim_one/extract_boundary_mesh_04.cc
@@ -0,0 +1,138 @@
+//----------------------------  extract_boundary_mesh_03.cc  ---------------------------
+//    $Id: bem.cc 22693 2010-11-11 20:11:47Z kanschat $
+//    Version: $Name$
+//
+//    Copyright (C) 2010 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  extract_boundary_mesh_03.cc  ---------------------------
+
+/*
+ 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.
+
+ Like _03, but set volume boundary indicator to 1 only on the upper
+ half of the sphere .
+*/
+
+#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>
+void 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)
+{
+  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;
+
+  for (; cell!=endc; ++cell){
+
+    face = surface_to_volume_mapping [cell];
+    Assert (face->at_boundary(), ExcInternalError());
+
+    deallog << "Surface cell: " << cell << " with vertices:" << std::endl;
+    for (unsigned int k=0; k<GeometryInfo<s_dim>::vertices_per_cell; ++k)
+      {
+       deallog << "  " << cell->vertex(k) << std::endl;
+       Assert (std::fabs(cell->vertex(k).distance (Point<spacedim>()) - 1) < 1e-12,
+               ExcInternalError());
+      }
+
+    deallog << "Volume face: " << face << " with vertices:" << std::endl;
+    for (unsigned int k=0; k<GeometryInfo<s_dim>::vertices_per_cell; ++k)
+      {
+       deallog << "  " << face->vertex(k) << std::endl;
+       Assert (std::fabs(face->vertex(k).distance (Point<spacedim>()) - 1) < 1e-12,
+               ExcInternalError());
+      }
+
+
+    for (unsigned int k=0; k<GeometryInfo<s_dim>::vertices_per_cell; ++k)
+      {
+       Point<spacedim> diff(face->vertex(k));
+       diff -= cell->vertex(k);
+       Assert (diff.square() == 0, ExcInternalError());
+      }
+  }
+}
+
+template <int dim, int spacedim>
+void save_mesh(const Triangulation<dim,spacedim>& tria)
+{
+  GridOut grid_out;
+  grid_out.write_gnuplot (tria, deallog.get_file_stream());
+}
+
+
+int main ()
+{
+
+  ofstream logfile("extract_boundary_mesh_04/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);
+    for (Triangulation<dim>::active_cell_iterator
+          cell = volume_mesh.begin_active();
+        cell != volume_mesh.end(); ++cell)
+      for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+       if (cell->at_boundary(f) && (cell->center()[dim-1]>0))
+         cell->face(f)->set_all_boundary_indicators (1);
+    
+    volume_mesh.set_boundary (1, boundary_description);
+    volume_mesh.set_boundary (0, boundary_description);
+    volume_mesh.refine_global (3);
+
+    const HyperBallBoundary<dim-1,dim> surface_description;
+    Triangulation<dim-1,dim> boundary_mesh;
+    boundary_mesh.set_boundary (1, surface_description);
+    boundary_mesh.set_boundary (0, surface_description);
+   
+    set<unsigned char> boundary_ids;
+    boundary_ids.insert(1);
+    
+    GridTools::extract_boundary_mesh (volume_mesh, boundary_mesh,
+                                     surface_to_volume_mapping,
+                                     boundary_ids);
+
+    deallog << volume_mesh.n_active_cells () << std::endl;
+    deallog << boundary_mesh.n_active_cells () << std::endl;
+    //save_mesh(boundary_mesh);
+
+
+    test_vertices_orientation(boundary_mesh, surface_to_volume_mapping);
+    //save_mesh(boundary_mesh);
+  }
+
+
+   return 0;
+ }