libMesh
sfc_partitioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // Local Includes
21 #include "libmesh/libmesh_config.h"
22 #include "libmesh/mesh_base.h"
23 #include "libmesh/sfc_partitioner.h"
24 #include "libmesh/libmesh_logging.h"
25 #include "libmesh/elem.h"
26 
27 #ifdef LIBMESH_HAVE_SFCURVES
28 namespace Sfc {
29 extern "C" {
30 # include "sfcurves.h"
31 }
32 }
33 #else
34 # include "libmesh/linear_partitioner.h"
35 #endif
36 
37 namespace libMesh
38 {
39 
40 
41 void SFCPartitioner::partition_range(MeshBase & mesh,
44  unsigned int n)
45 {
46  libmesh_assert_greater (n, 0);
47 
48  // Check for an easy return
49  if (n == 1)
50  {
51  this->single_partition_range (beg, end);
52  return;
53  }
54 
55  // What to do if the sfcurves library IS NOT present
56 #ifndef LIBMESH_HAVE_SFCURVES
57 
58  libmesh_here();
59  libMesh::err << "ERROR: The library has been built without" << std::endl
60  << "Space Filling Curve support. Using a linear" << std::endl
61  << "partitioner instead!" << std::endl;
62 
64  lp.partition_range (mesh, beg, end, n);
65 
66  // What to do if the sfcurves library IS present
67 #else
68 
69  LOG_SCOPE("partition_range()", "SFCPartitioner");
70 
71  const dof_id_type n_range_elem = std::distance(beg, end);
72  const dof_id_type n_elem = mesh.n_elem();
73 
74  // The forward_map maps the range's element ids into a contiguous
75  // block of indices.
76  std::vector<dof_id_type> forward_map (n_elem, DofObject::invalid_id);
77 
78  // the reverse_map maps the contiguous ids back
79  // to active elements
80  std::vector<Elem *> reverse_map (n_range_elem, libmesh_nullptr);
81 
82  std::vector<double> x (n_range_elem);
83  std::vector<double> y (n_range_elem);
84  std::vector<double> z (n_range_elem);
85  std::vector<int> table (n_range_elem);
86 
87  // Map the range's element ids into a contiguous range.
88  {
90  dof_id_type el_num = 0;
91 
92  for (; it != end; ++it)
93  {
94  libmesh_assert_less ((*it)->id(), forward_map.size());
95  libmesh_assert_less (el_num, reverse_map.size());
96 
97  forward_map[(*it)->id()] = el_num;
98  reverse_map[el_num] = *it;
99  el_num++;
100  }
101  libmesh_assert_equal_to (el_num, n_range_elem);
102  }
103 
104  // Get the centroid for each range element.
105  {
107 
108  for (; it != end; ++it)
109  {
110  const Elem * elem = *it;
111 
112  libmesh_assert_less (elem->id(), forward_map.size());
113 
114  const Point p = elem->centroid();
115 
116  x[forward_map[elem->id()]] = p(0);
117  y[forward_map[elem->id()]] = p(1);
118  z[forward_map[elem->id()]] = p(2);
119  }
120  }
121 
122  // We need an integer reference to pass to the Sfc interface.
123  int size = static_cast<int>(n_range_elem);
124 
125  // build the space-filling curve
126  if (_sfc_type == "Hilbert")
127  Sfc::hilbert (&x[0], &y[0], &z[0], &size, &table[0]);
128 
129  else if (_sfc_type == "Morton")
130  Sfc::morton (&x[0], &y[0], &z[0], &size, &table[0]);
131 
132  else
133  {
134  libmesh_here();
135  libMesh::err << "ERROR: Unknown type: " << _sfc_type << std::endl
136  << " Valid types are" << std::endl
137  << " \"Hilbert\"" << std::endl
138  << " \"Morton\"" << std::endl
139  << " " << std::endl
140  << "Proceeding with a Hilbert curve." << std::endl;
141 
142  Sfc::hilbert (&x[0], &y[0], &z[0], &size, &table[0]);
143  }
144 
145 
146  // Assign the partitioning to the range elements
147  {
148  // {
149  // std::ofstream out ("sfc.dat");
150  // out << "variables=x,y,z" << std::endl;
151  // out << "zone f=point" << std::endl;
152  // for (unsigned int i=0; i<n_range_elem; i++)
153  // out << x[i] << " " << y[i] << " " << z[i] << std::endl;
154  // }
155 
156  const dof_id_type blksize = (n_range_elem + n - 1) / n;
157 
158  for (dof_id_type i=0; i<n_range_elem; i++)
159  {
160  libmesh_assert_less (static_cast<unsigned int>(table[i] - 1), reverse_map.size());
161 
162  Elem * elem = reverse_map[table[i] - 1];
163 
164  elem->processor_id() = cast_int<processor_id_type>(i/blksize);
165  }
166  }
167 
168 #endif
169 }
170 
171 
172 
173 void SFCPartitioner::_do_partition (MeshBase & mesh,
174  const unsigned int n)
175 {
176  this->partition_range(mesh,
177  mesh.active_elements_begin(),
178  mesh.active_elements_end(),
179  n);
180 }
181 
182 } // namespace libMesh
The definition of the element_iterator struct.
Definition: mesh_base.h:1476
OStreamProxy err
dof_id_type n_elem(const MeshBase::const_element_iterator &begin, const MeshBase::const_element_iterator &end)
Count up the number of elements of a specific type (as defined by an iterator range).
Definition: mesh_tools.C:656
The LinearPartitioner simply takes the element list and splits it into equal-sized chunks assigned to...
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
MeshBase & mesh
const class libmesh_nullptr_t libmesh_nullptr
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
The libMesh namespace provides an interface to certain functionality in the library.
virtual void partition_range(MeshBase &mesh, MeshBase::element_iterator it, MeshBase::element_iterator end, const unsigned int n) libmesh_override
Called by the SubdomainPartitioner to partition elements in the range (it, end).
Real distance(const Point &p)
This is the MeshBase class.
Definition: mesh_base.h:68
PetscErrorCode Vec x
virtual element_iterator active_elements_begin()=0
Active, local, and negation forms of the element iterators described above.
virtual element_iterator active_elements_end()=0
virtual Point centroid() const
Definition: elem.C:446
dof_id_type id() const
Definition: dof_object.h:632
virtual dof_id_type n_elem() const =0
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
uint8_t dof_id_type
Definition: id_types.h:64
processor_id_type processor_id() const
Definition: dof_object.h:694