From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 29 Oct 2001 08:43:35 +0000 (+0000)
Subject: Provide subdomain ids to triangulation cells.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6b55e5f182b90cadc94403c205a6a4e606082539;p=dealii-svn.git

Provide subdomain ids to triangulation cells.


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

diff --git a/deal.II/deal.II/include/grid/tria_accessor.h b/deal.II/deal.II/include/grid/tria_accessor.h
index 5385ceff5e..2027172307 100644
--- a/deal.II/deal.II/include/grid/tria_accessor.h
+++ b/deal.II/deal.II/include/grid/tria_accessor.h
@@ -14,6 +14,7 @@
 #define __deal2__tria_accessor_h
 
 
+#include <base/config.h>
 #include <base/exceptions.h>
 #include <grid/tria_iterator_base.h>
 
@@ -1966,14 +1967,28 @@ class CellAccessor :  public TriaObjectAccessor<dim,dim>
     face (const unsigned int i) const;
     
 				     /**
-				      * Return the material id of this cell.
+				      * Return the material id of this
+				      * cell.
 				      */
     unsigned char material_id () const;
 
 				     /**
-				      * Set the material id of this cell.
+				      * Set the material id of this
+				      * cell.
 				      */
-    void set_material_id (const unsigned char) const;
+    void set_material_id (const unsigned char new_material_id) const;
+
+				     /**
+				      * Return the subdomain id of
+				      * this cell.
+				      */
+    unsigned int subdomain_id () const;
+
+				     /**
+				      * Set the subdomain id of this
+				      * cell.
+				      */
+    void set_subdomain_id (const unsigned int new_subdomain_id) const;
 
 				     /**
 				      * Test whether the cell has children
@@ -2056,6 +2071,7 @@ template <> bool CellAccessor<3>::at_boundary () const;
 template <> unsigned char CellAccessor<3>::material_id () const;
 template <> void CellAccessor<3>::set_material_id (const unsigned char mat_id) const;
 template <> bool CellAccessor<3>::point_inside (const Point<3> &) const;
+
 template <> bool CellAccessor<1>::has_boundary_lines () const;
 
 template <> double TriaObjectAccessor<2, 2>::measure () const;
diff --git a/deal.II/deal.II/include/grid/tria_accessor.templates.h b/deal.II/deal.II/include/grid/tria_accessor.templates.h
index 367fdb0599..1fd585ed08 100644
--- a/deal.II/deal.II/include/grid/tria_accessor.templates.h
+++ b/deal.II/deal.II/include/grid/tria_accessor.templates.h
@@ -10,6 +10,7 @@
 #define __deal2__tria_accessor_templates_h
 
 
+#include <base/config.h>
 #include <grid/tria.h>
 #include <grid/tria_levels.h>
 #include <grid/tria_iterator.h>
diff --git a/deal.II/deal.II/include/grid/tria_levels.h b/deal.II/deal.II/include/grid/tria_levels.h
index 4c3c454c47..8e91ee985b 100644
--- a/deal.II/deal.II/include/grid/tria_levels.h
+++ b/deal.II/deal.II/include/grid/tria_levels.h
@@ -14,6 +14,7 @@
 #define __deal2__tria_levels_h
 
 
+#include <base/config.h>
 #include <vector>
 #include <grid/tria_line.h>
 #include <grid/tria_quad.h>
@@ -127,6 +128,18 @@ class TriangulationLevel<0>
 				      *  to the mother cell of this cell).
 				      */
     std::vector<std::pair<int,int> > neighbors;
+
+				     /**
+				      * One integer per cell to store
+				      * which subdomain it belongs
+				      * to. This field is most often
+				      * used in parallel computations,
+				      * where it denotes which
+				      * processor shall work on the
+				      * cells with a given subdomain
+				      * number.
+				      */
+    std::vector<unsigned int> subdomain_ids;
     
 				     /**
 				      *  Reserve enough space to accomodate
diff --git a/deal.II/deal.II/source/grid/tria_accessor.cc b/deal.II/deal.II/source/grid/tria_accessor.cc
index c1c548281d..9701c1dd08 100644
--- a/deal.II/deal.II/source/grid/tria_accessor.cc
+++ b/deal.II/deal.II/source/grid/tria_accessor.cc
@@ -1562,7 +1562,7 @@ void CellAccessor<1>::set_material_id (const unsigned char mat_id) const
 {
   Assert (used(), TriaAccessor<1>::ExcCellNotUsed());
   tria->levels[present_level]->lines.material_id[present_index]
-    = mat_id;						 
+    = mat_id;
 };
 
 
@@ -1705,6 +1705,26 @@ bool CellAccessor<3>::point_inside (const Point<3> &) const
 /*------------------------ Functions: CellAccessor<dim> -----------------------*/
 
 
+template <int dim>
+unsigned int CellAccessor<dim>::subdomain_id () const
+{
+  Assert (used(), typename TriaAccessor<dim>::ExcCellNotUsed());
+  return tria->levels[present_level]->subdomain_ids[present_index];
+};
+
+
+
+template <int dim>
+void
+CellAccessor<dim>::set_subdomain_id (const unsigned int new_subdomain_id) const
+{
+  Assert (used(), typename TriaAccessor<dim>::ExcCellNotUsed());
+  tria->levels[present_level]->subdomain_ids[present_index]
+    = new_subdomain_id;
+};
+
+
+
 template <int dim>
 void CellAccessor<dim>::set_neighbor (const unsigned int i,
 				      const TriaIterator<dim,CellAccessor<dim> > &pointer) const
diff --git a/tests/deal.II/Makefile b/tests/deal.II/Makefile
index f421960a69..27474438a4 100644
--- a/tests/deal.II/Makefile
+++ b/tests/deal.II/Makefile
@@ -45,12 +45,14 @@ filtered_matrix.exe          : filtered_matrix.go          $(lib-1d) $(lib-2d) $
 boundaries.exe               : boundaries.go               $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 sparsity_pattern.exe         : sparsity_pattern.go         $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 grid_tools.exe               : grid_tools.go               $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
+subdomain_ids.exe            : subdomain_ids.go            $(lib-1d) $(lib-2d) $(lib-3d) $(libraries)
 
 tests = grid_test grid_transform dof_test data_out derivatives gradients \
 	constraints block_matrices second_derivatives derivative_approximation \
         matrices error_estimator intergrid_constraints intergrid_map \
         wave-test-3 dof_renumbering support_point_map filtered_matrix \
-	boundaries sparsity_pattern grid_tools
+	boundaries sparsity_pattern grid_tools subdomain_ids
+
 ############################################################
 
 
diff --git a/tests/deal.II/subdomain_ids.cc b/tests/deal.II/subdomain_ids.cc
new file mode 100644
index 0000000000..be7d3ecc87
--- /dev/null
+++ b/tests/deal.II/subdomain_ids.cc
@@ -0,0 +1,185 @@
+//----------------------------  subdomain_ids.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 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.
+//
+//----------------------------  subdomain_ids.cc  ---------------------------
+
+// check some things about subdomain_ids
+
+
+#include <base/logstream.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_iterator.h>
+#include <dofs/dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_dgq.h>
+#include <fe/fe_q.h>
+#include <dofs/dof_tools.h>
+
+#include <fstream>
+#include <algorithm>
+#include <cmath>
+
+
+std::ofstream logfile("subdomain_ids.output");
+
+
+DeclException2 (ExcNumberMismatch,
+		int, int,
+		<< "The numbers " << arg1 << " and " << arg2
+		<< " should be equation, but are not.");
+
+
+template <int dim>
+void test ()
+{
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria, -1, 1);
+  tria.refine_global (3);
+
+				   // we now have a number of cells,
+				   // flag them with some subdomain
+				   // ids based on their position, in
+				   // particular we take the quadrant
+				   // (octant)
+  typename Triangulation<dim>::active_cell_iterator
+    cell = tria.begin_active (),
+    endc = tria.end ();
+  for (; cell!=endc; ++cell)
+    {
+      unsigned int subdomain = 0;
+      for (unsigned int d=0; d<dim; ++d)
+	if (cell->center()(d) > 0)
+	  subdomain |= (1<<d);
+      Assert (subdomain < (1<<dim), ExcInternalError());
+
+      cell->set_subdomain_id (subdomain);
+    };
+
+				   // refine twice (for 3d only once)
+				   // and count again the numbers of
+				   // cells in each subdomain
+  tria.refine_global ((dim != 3) ? 2 : 1);
+  if (true)
+    {
+      cell = tria.begin_active();
+      endc = tria.end();
+      std::vector<unsigned int> subdomain_cells(1<<dim, 0);
+      for (; cell!=endc; ++cell)
+	{
+	  Assert (cell->subdomain_id() < (1<<dim), ExcInternalError());
+	  ++subdomain_cells[cell->subdomain_id()];
+	};
+      for (unsigned int i=0; i<(1<<dim); ++i)
+	Assert (subdomain_cells[i] == tria.n_active_cells()/(1<<dim),
+		ExcNumberMismatch(subdomain_cells[i],
+				  tria.n_active_cells()/(1<<dim)));
+      deallog << "Check 1 (dim=" << dim << ") ok" << std::endl;
+    };
+  
+				   // coarsen once and check again
+  if (true)
+    {
+      cell = tria.begin_active();
+      endc = tria.end();
+      for (; cell!=endc; ++cell)
+	cell->set_coarsen_flag ();
+      tria.execute_coarsening_and_refinement ();
+      
+      cell = tria.begin_active();
+      endc = tria.end();
+      std::vector<unsigned int> subdomain_cells(1<<dim, 0);
+      for (; cell!=endc; ++cell)
+	{
+	  Assert (cell->subdomain_id() < (1<<dim), ExcInternalError());
+	  ++subdomain_cells[cell->subdomain_id()];
+	};
+      for (unsigned int i=0; i<(1<<dim); ++i)
+	Assert (subdomain_cells[i] == tria.n_active_cells()/(1<<dim),
+		ExcNumberMismatch(subdomain_cells[i],
+				  tria.n_active_cells()/(1<<dim)));
+      deallog << "Check 2 (dim=" << dim << ") ok" << std::endl;
+    };
+
+				   // check 3: assign DQ2 elements to
+				   // cells and count degrees of
+				   // freedom associated with each
+				   // subdomain
+  if (true)
+    {
+      FE_DGQ<dim> fe(2);
+      DoFHandler<dim> dof_handler (tria);
+      dof_handler.distribute_dofs (fe);
+      std::vector<bool> selected_dofs (dof_handler.n_dofs());
+      for (unsigned int subdomain=0; subdomain<(1<<dim); ++subdomain)
+	{
+	  DoFTools::extract_subdomain_dofs (dof_handler, subdomain,
+					    selected_dofs);
+	  Assert (static_cast<unsigned int>(std::count (selected_dofs.begin(),
+							selected_dofs.end(),
+							true))
+		  ==
+		  dof_handler.n_dofs() / (1<<dim),
+		  ExcNumberMismatch(std::count (selected_dofs.begin(),
+						selected_dofs.end(),
+						true),
+				    dof_handler.n_dofs() / (1<<dim)));
+	}
+      deallog << "Check 3 (dim=" << dim << ") ok" << std::endl;      
+    };
+
+
+				   // check 4: check again for
+				   // continuous elements. note that
+				   // the number of dofs here is
+				   // different, since dofs can be on
+				   // several subdomain at once
+  if (true)
+    {
+      FE_Q<dim> fe(1);
+      DoFHandler<dim> dof_handler (tria);
+      dof_handler.distribute_dofs (fe);
+
+      const unsigned int cells_per_direction
+	= static_cast<unsigned int>(rint(pow(tria.n_active_cells(), 1./dim)));
+      
+      std::vector<bool> selected_dofs (dof_handler.n_dofs());
+      for (unsigned int subdomain=0; subdomain<(1<<dim); ++subdomain)
+	{
+	  DoFTools::extract_subdomain_dofs (dof_handler, subdomain,
+					    selected_dofs);
+	  Assert (static_cast<unsigned int>(std::count (selected_dofs.begin(),
+							selected_dofs.end(),
+							true))
+		  ==
+		  pow(cells_per_direction/2+1,dim),
+		  ExcNumberMismatch(std::count (selected_dofs.begin(),
+						selected_dofs.end(),
+						true),
+				    pow(cells_per_direction/2+1,dim)));
+	}
+      deallog << "Check 4 (dim=" << dim << ") ok" << std::endl;      
+    };
+};
+
+
+int main ()
+{
+  logfile.precision(4);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+
+  test<1> ();
+  test<2> ();
+  test<3> ();
+  
+  return 0;
+};
diff --git a/tests/deal.II/subdomain_ids.checked b/tests/deal.II/subdomain_ids.checked
new file mode 100644
index 0000000000..b8cc67cb10
--- /dev/null
+++ b/tests/deal.II/subdomain_ids.checked
@@ -0,0 +1,13 @@
+
+DEAL::Check 1 (dim=1) ok
+DEAL::Check 2 (dim=1) ok
+DEAL::Check 3 (dim=1) ok
+DEAL::Check 4 (dim=1) ok
+DEAL::Check 1 (dim=2) ok
+DEAL::Check 2 (dim=2) ok
+DEAL::Check 3 (dim=2) ok
+DEAL::Check 4 (dim=2) ok
+DEAL::Check 1 (dim=3) ok
+DEAL::Check 2 (dim=3) ok
+DEAL::Check 3 (dim=3) ok
+DEAL::Check 4 (dim=3) ok