Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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/elem.h"
23 : #include "libmesh/enum_partitioner_type.h"
24 : #include "libmesh/libmesh_logging.h"
25 : #include "libmesh/mesh_base.h"
26 : #include "libmesh/sfc_partitioner.h"
27 :
28 : #ifdef LIBMESH_HAVE_SFCURVES
29 : namespace Sfc {
30 : extern "C" {
31 : # include "sfcurves.h"
32 : }
33 : }
34 : #else
35 : # include "libmesh/linear_partitioner.h"
36 : #endif
37 :
38 : namespace libMesh
39 : {
40 :
41 :
42 142 : PartitionerType SFCPartitioner::type() const
43 : {
44 142 : return SFC_PARTITIONER;
45 : }
46 :
47 :
48 813 : void SFCPartitioner::partition_range(MeshBase & mesh,
49 : MeshBase::element_iterator beg,
50 : MeshBase::element_iterator end,
51 : unsigned int n)
52 : {
53 : // Check for easy returns
54 813 : if (beg == end)
55 213 : return;
56 :
57 600 : if (n == 1)
58 : {
59 0 : this->single_partition_range (beg, end);
60 0 : return;
61 : }
62 :
63 102 : libmesh_assert_greater (n, 0);
64 :
65 : // What to do if the sfcurves library IS NOT present
66 : #ifndef LIBMESH_HAVE_SFCURVES
67 :
68 396 : libmesh_do_once(
69 : libMesh::out << "ERROR: The library has been built without" << std::endl
70 : << "Space Filling Curve support. Using a linear" << std::endl
71 : << "partitioner instead!" << std::endl;);
72 :
73 : LinearPartitioner lp;
74 792 : lp.partition_range (mesh, beg, end, n);
75 :
76 : // What to do if the sfcurves library IS present
77 : #else
78 :
79 204 : LOG_SCOPE("partition_range()", "SFCPartitioner");
80 :
81 : // We don't yet support distributed meshes with this Partitioner
82 204 : if (!mesh.is_serial())
83 0 : libmesh_not_implemented();
84 :
85 204 : const dof_id_type n_range_elem = std::distance(beg, end);
86 204 : const dof_id_type max_elem_id = mesh.max_elem_id();
87 :
88 : // The forward_map maps the range's element ids into a contiguous
89 : // block of indices.
90 306 : std::vector<dof_id_type> forward_map (max_elem_id, DofObject::invalid_id);
91 :
92 : // the reverse_map maps the contiguous ids back
93 : // to active elements
94 306 : std::vector<Elem *> reverse_map (n_range_elem, nullptr);
95 :
96 306 : std::vector<double> x (n_range_elem);
97 306 : std::vector<double> y (n_range_elem);
98 306 : std::vector<double> z (n_range_elem);
99 306 : std::vector<int> table (n_range_elem);
100 :
101 : // Map the range's element ids into a contiguous range.
102 102 : dof_id_type el_num = 0;
103 :
104 11070 : for (auto & elem : as_range(beg, end))
105 : {
106 5382 : libmesh_assert_less (elem->id(), forward_map.size());
107 5382 : libmesh_assert_less (el_num, reverse_map.size());
108 :
109 10764 : forward_map[elem->id()] = el_num;
110 10764 : reverse_map[el_num] = elem;
111 10764 : el_num++;
112 : }
113 102 : libmesh_assert_equal_to (el_num, n_range_elem);
114 :
115 : // Get the vertex average for each range element.
116 11070 : for (const auto & elem : as_range(beg, end))
117 : {
118 5382 : libmesh_assert_less (elem->id(), forward_map.size());
119 :
120 10764 : const Point p = elem->vertex_average();
121 :
122 16146 : x[forward_map[elem->id()]] = double(p(0));
123 10764 : y[forward_map[elem->id()]] = double(p(1));
124 16146 : z[forward_map[elem->id()]] = double(p(2));
125 : }
126 :
127 : // We need an integer reference to pass to the Sfc interface.
128 204 : int size = static_cast<int>(n_range_elem);
129 :
130 : // build the space-filling curve
131 204 : if (_sfc_type == "Hilbert")
132 136 : Sfc::hilbert (x.data(),
133 : y.data(),
134 : z.data(),
135 : &size,
136 : table.data());
137 :
138 68 : else if (_sfc_type == "Morton")
139 68 : Sfc::morton (x.data(),
140 : y.data(),
141 : z.data(),
142 : &size,
143 : table.data());
144 :
145 : else
146 : {
147 0 : libMesh::out << "ERROR: Unknown type: " << _sfc_type << std::endl
148 0 : << " Valid types are" << std::endl
149 0 : << " \"Hilbert\"" << std::endl
150 0 : << " \"Morton\"" << std::endl
151 0 : << " " << std::endl
152 0 : << "Partitioning with a Hilbert curve." << std::endl;
153 :
154 0 : Sfc::hilbert (x.data(),
155 : y.data(),
156 : z.data(),
157 : &size,
158 : table.data());
159 : }
160 :
161 :
162 : // Assign the partitioning to the range elements
163 : {
164 : // {
165 : // std::ofstream out ("sfc.dat");
166 : // out << "variables=x,y,z" << std::endl;
167 : // out << "zone f=point" << std::endl;
168 : // for (unsigned int i=0; i<n_range_elem; i++)
169 : // out << x[i] << " " << y[i] << " " << z[i] << std::endl;
170 : // }
171 :
172 204 : const dof_id_type blksize = (n_range_elem + n - 1) / n;
173 :
174 10968 : for (dof_id_type i=0; i<n_range_elem; i++)
175 : {
176 5382 : libmesh_assert_less (table[i] - 1, reverse_map.size());
177 :
178 16146 : Elem * elem = reverse_map[table[i] - 1];
179 :
180 10764 : elem->processor_id() = cast_int<processor_id_type>(i/blksize);
181 : }
182 : }
183 :
184 : #endif
185 : }
186 :
187 :
188 :
189 813 : void SFCPartitioner::_do_partition (MeshBase & mesh,
190 : const unsigned int n)
191 : {
192 813 : this->partition_range(mesh,
193 1626 : mesh.active_elements_begin(),
194 921 : mesh.active_elements_end(),
195 216 : n);
196 813 : }
197 :
198 : } // namespace libMesh
|
