sfc_partitioner.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 // C++ Includes -----------------------------------
21 
22 // Local Includes -----------------------------------
23 #include "libmesh/libmesh_config.h"
24 #include "libmesh/mesh_base.h"
27 #include "libmesh/elem.h"
28 
29 #ifdef LIBMESH_HAVE_SFCURVES
30  namespace Sfc {
31  extern "C" {
32 # include "sfcurves.h"
33  }
34  }
35 #else
37 #endif
38 
39 namespace libMesh
40 {
41 
42 
43 // ------------------------------------------------------------
44 // SFCPartitioner implementation
46  const unsigned int n)
47 {
48 
50 
51  // Check for an easy return
52  if (n == 1)
53  {
54  this->single_partition (mesh);
55  return;
56  }
57 
58 // What to do if the sfcurves library IS NOT present
59 #ifndef LIBMESH_HAVE_SFCURVES
60 
61  libmesh_here();
62  libMesh::err << "ERROR: The library has been built without" << std::endl
63  << "Space Filling Curve support. Using a linear" << std::endl
64  << "partitioner instead!" << std::endl;
65 
67 
68  lp.partition (mesh, n);
69 
70 // What to do if the sfcurves library IS present
71 #else
72 
73  START_LOG("sfc_partition()", "SFCPartitioner");
74 
75  const dof_id_type n_active_elem = mesh.n_active_elem();
76  const dof_id_type n_elem = mesh.n_elem();
77 
78  // the forward_map maps the active element id
79  // into a contiguous block of indices
80  std::vector<dof_id_type>
81  forward_map (n_elem, DofObject::invalid_id);
82 
83  // the reverse_map maps the contiguous ids back
84  // to active elements
85  std::vector<Elem*> reverse_map (n_active_elem, NULL);
86 
87  int size = static_cast<int>(n_active_elem);
88  std::vector<double> x (size);
89  std::vector<double> y (size);
90  std::vector<double> z (size);
91  std::vector<int> table (size);
92 
93 
94  // We need to map the active element ids into a
95  // contiguous range.
96  {
97 // active_elem_iterator elem_it (mesh.elements_begin());
98 // const active_elem_iterator elem_end(mesh.elements_end());
99 
101  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
102 
103  dof_id_type el_num = 0;
104 
105  for (; elem_it != elem_end; ++elem_it)
106  {
107  libmesh_assert_less ((*elem_it)->id(), forward_map.size());
108  libmesh_assert_less (el_num, reverse_map.size());
109 
110  forward_map[(*elem_it)->id()] = el_num;
111  reverse_map[el_num] = *elem_it;
112  el_num++;
113  }
114  libmesh_assert_equal_to (el_num, n_active_elem);
115  }
116 
117 
118  // Get the centroid for each active element
119  {
120 // const_active_elem_iterator elem_it (mesh.const_elements_begin());
121 // const const_active_elem_iterator elem_end(mesh.const_elements_end());
122 
124  const MeshBase::element_iterator elem_end = mesh.active_elements_end();
125 
126  for (; elem_it != elem_end; ++elem_it)
127  {
128  const Elem* elem = *elem_it;
129 
130  libmesh_assert_less (elem->id(), forward_map.size());
131 
132  const Point p = elem->centroid();
133 
134  x[forward_map[elem->id()]] = p(0);
135  y[forward_map[elem->id()]] = p(1);
136  z[forward_map[elem->id()]] = p(2);
137  }
138  }
139 
140  // build the space-filling curve
141  if (_sfc_type == "Hilbert")
142  Sfc::hilbert (&x[0], &y[0], &z[0], &size, &table[0]);
143 
144  else if (_sfc_type == "Morton")
145  Sfc::morton (&x[0], &y[0], &z[0], &size, &table[0]);
146 
147  else
148  {
149  libmesh_here();
150  libMesh::err << "ERROR: Unknown type: " << _sfc_type << std::endl
151  << " Valid types are" << std::endl
152  << " \"Hilbert\"" << std::endl
153  << " \"Morton\"" << std::endl
154  << " " << std::endl
155  << "Proceeding with a Hilbert curve." << std::endl;
156 
157  Sfc::hilbert (&x[0], &y[0], &z[0], &size, &table[0]);
158  }
159 
160 
161  // Assign the partitioning to the active elements
162  {
163 // {
164 // std::ofstream out ("sfc.dat");
165 // out << "variables=x,y,z" << std::endl;
166 // out << "zone f=point" << std::endl;
167 
168 // for (unsigned int i=0; i<n_active_elem; i++)
169 // out << x[i] << " "
170 // << y[i] << " "
171 // << z[i] << std::endl;
172 // }
173 
174  const dof_id_type blksize = (n_active_elem+n-1)/n;
175 
176  for (dof_id_type i=0; i<n_active_elem; i++)
177  {
178  libmesh_assert_less (static_cast<unsigned int>(table[i]-1), reverse_map.size());
179 
180  Elem* elem = reverse_map[table[i]-1];
181 
182  elem->processor_id() = libmesh_cast_int<processor_id_type>
183  (i/blksize);
184  }
185  }
186 
187  STOP_LOG("sfc_partition()", "SFCPartitioner");
188 
189 #endif
190 
191 }
192 
193 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo