Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "KokkosMesh.h"
11 :
12 : #include "Assembly.h"
13 : #include "MooseMesh.h"
14 :
15 : #include "libmesh/elem_side_builder.h"
16 : #include "libmesh/reference_elem.h"
17 :
18 : namespace Moose::Kokkos
19 : {
20 :
21 : void
22 2460 : Mesh::update()
23 : {
24 2460 : initMap();
25 2460 : initElement();
26 :
27 2460 : _initialized = true;
28 :
29 2460 : if (_needs_element_geometry)
30 181 : initElementGeometry();
31 2460 : if (_needs_element_side_geometry)
32 181 : initElementSideGeometry();
33 2460 : }
34 :
35 : void
36 2460 : Mesh::initMap()
37 : {
38 2460 : if (!_maps)
39 2460 : _maps = std::make_shared<MeshMap>();
40 :
41 2460 : if (_elem_id_integer == libMesh::invalid_uint)
42 7380 : _elem_id_integer = _mesh.getMesh().add_elem_integer("kokkos_contiguous_elem_id");
43 :
44 2460 : if (_node_id_integer == libMesh::invalid_uint)
45 7380 : _node_id_integer = _mesh.getMesh().add_node_integer("kokkos_contiguous_node_id");
46 :
47 2460 : _num_local_elems = 0;
48 2460 : _num_ghost_elems = 0;
49 2460 : _num_local_nodes = 0;
50 :
51 2460 : _maps->subdomain_id_mapping.clear();
52 2460 : _maps->boundary_id_mapping.clear();
53 2460 : _maps->elem_type_id_mapping.clear();
54 2460 : _maps->local_nodes.clear();
55 2460 : _maps->ghost_node_id_mapping.clear();
56 2460 : _maps->ghost_elem_id_mapping.clear();
57 2460 : _maps->subdomain_elem_id_ranges.clear();
58 2460 : _maps->subdomain_node_ids.clear();
59 2460 : _maps->boundary_node_ids.clear();
60 :
61 2460 : std::unordered_set<ElemType> elem_types;
62 2460 : std::unordered_set<Node *> ghost_nodes;
63 :
64 5621 : for (const auto subdomain : _mesh.meshSubdomains())
65 : {
66 3161 : _maps->subdomain_id_mapping[subdomain] = _maps->subdomain_id_mapping.size();
67 3161 : _maps->subdomain_node_ids[subdomain];
68 :
69 3161 : dof_id_type begin = libMesh::DofObject::invalid_id;
70 3161 : dof_id_type eid = libMesh::DofObject::invalid_id;
71 :
72 668509 : for (auto elem : _mesh.getMesh().active_local_subdomain_elements_ptr_range(subdomain))
73 : {
74 665348 : elem_types.insert(elem->type());
75 :
76 665348 : eid = _num_local_elems;
77 665348 : elem->set_extra_integer(_elem_id_integer, _num_local_elems++);
78 :
79 665348 : if (begin == libMesh::DofObject::invalid_id)
80 2969 : begin = eid;
81 :
82 3357034 : for (auto & node : elem->node_ref_range())
83 2691686 : if (node.processor_id() != _mesh.processor_id())
84 13938 : ghost_nodes.insert(&node);
85 3161 : }
86 :
87 3161 : _maps->subdomain_elem_id_ranges[subdomain] = eid != libMesh::DofObject::invalid_id
88 3418 : ? std::make_pair(begin, eid + 1)
89 1805 : : std::make_pair(begin, eid);
90 : }
91 :
92 12117 : for (const auto boundary : _mesh.meshBoundaryIds())
93 : {
94 9657 : _maps->boundary_id_mapping[boundary] = _maps->boundary_id_mapping.size();
95 9657 : _maps->boundary_node_ids[boundary];
96 : }
97 :
98 : // Enumerate one off-process neighbor layer for element-side geometry data.
99 2460 : if (_needs_element_side_geometry)
100 21554 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
101 102487 : for (unsigned int side = 0; side < elem->n_sides(); ++side)
102 : {
103 81114 : const auto neighbor = elem->neighbor_ptr(side);
104 79330 : if (neighbor && neighbor != libMesh::remote_elem &&
105 79408 : neighbor->processor_id() != _mesh.processor_id() &&
106 16 : !_maps->ghost_elem_id_mapping.count(neighbor))
107 : {
108 32 : _maps->ghost_elem_id_mapping[neighbor] = _num_local_elems + _num_ghost_elems++;
109 32 : elem_types.insert(neighbor->type());
110 : }
111 181 : }
112 :
113 : // Setup node maps
114 :
115 4916 : for (const auto type : elem_types)
116 2456 : _maps->elem_type_id_mapping[type] = _maps->elem_type_id_mapping.size();
117 :
118 3790 : _maps->local_nodes.insert(
119 2660 : _maps->local_nodes.end(), _mesh.localNodesBegin(), _mesh.localNodesEnd());
120 2460 : _maps->local_nodes.insert(_maps->local_nodes.end(), ghost_nodes.begin(), ghost_nodes.end());
121 :
122 734666 : for (auto node : _maps->local_nodes)
123 : {
124 732206 : if (node->processor_id() == _mesh.processor_id())
125 : {
126 724774 : node->set_extra_integer(_node_id_integer, _num_local_nodes);
127 :
128 1535881 : for (auto subdomain : _mesh.getNodeBlockIds(*node))
129 811107 : _maps->subdomain_node_ids[subdomain].push_back(_num_local_nodes);
130 : }
131 : else
132 7432 : _maps->ghost_node_id_mapping[node] = _num_local_nodes;
133 :
134 732206 : ++_num_local_nodes;
135 : }
136 :
137 126507 : for (const auto bnd_node : as_range(_mesh.bndNodesBegin(), _mesh.bndNodesEnd()))
138 124047 : if (bnd_node->_node->processor_id() == _mesh.processor_id())
139 105817 : _maps->boundary_node_ids[bnd_node->_bnd_id].insert(getContiguousNodeID(bnd_node->_node));
140 2460 : }
141 :
142 : void
143 2460 : Mesh::initElement()
144 : {
145 : // Cache reference element data
146 :
147 2460 : const auto num_elem_types = getNumLocalElementTypes();
148 :
149 2460 : _num_sides.create(num_elem_types);
150 2460 : _num_nodes.create(num_elem_types);
151 2460 : _num_side_nodes.create(num_elem_types);
152 2460 : _local_side_node.create(num_elem_types);
153 :
154 4916 : for (auto & [elem_type, elem_type_id] : getElementTypeMap())
155 : {
156 2456 : auto elem = &libMesh::ReferenceElem::get(elem_type);
157 :
158 2456 : _num_sides[elem_type_id] = elem->n_sides();
159 2456 : _num_nodes[elem_type_id] = elem->n_nodes();
160 :
161 2456 : _num_side_nodes[elem_type_id].create(elem->n_sides());
162 2456 : _local_side_node[elem_type_id].create(elem->n_nodes(), elem->n_sides());
163 :
164 11740 : for (unsigned int side = 0; side < elem->n_sides(); ++side)
165 : {
166 9284 : _num_side_nodes[elem_type_id][side] = elem->side_ptr(side)->n_nodes();
167 :
168 28162 : for (unsigned int node = 0; node < elem->side_ptr(side)->n_nodes(); ++node)
169 18878 : _local_side_node[elem_type_id](node, side) = elem->local_side_node(side, node);
170 : }
171 :
172 2456 : _num_side_nodes[elem_type_id].moveToDevice();
173 2456 : _local_side_node[elem_type_id].moveToDevice();
174 : }
175 :
176 2460 : _num_sides.moveToDevice();
177 2460 : _num_nodes.moveToDevice();
178 2460 : _num_side_nodes.copyToDevice();
179 2460 : _local_side_node.copyToDevice();
180 :
181 : // Cache element data
182 :
183 1130 : const auto num_local_plus_one_neighbor_layer_elems =
184 1330 : getNumLocalAndPossiblyOneNeighborLayerGhostElements();
185 2460 : const auto num_elems = getNumLocalElements();
186 2460 : const auto num_subdomains = getNumSubdomains();
187 :
188 2460 : _elem_info.create(num_local_plus_one_neighbor_layer_elems);
189 2460 : _elem_neighbor.create(_mesh.getMaxSidesPerElem(), num_elems);
190 2460 : _extra_elem_ids.create(num_elems, _mesh.getMesh().n_elem_integers());
191 2460 : _starting_elem_id.create(num_subdomains);
192 :
193 2460 : _elem_neighbor = libMesh::DofObject::invalid_id;
194 2460 : _extra_elem_ids = libMesh::DofObject::invalid_id;
195 :
196 5621 : for (const auto subdomain : _mesh.meshSubdomains())
197 : {
198 3161 : const auto sid = getContiguousSubdomainID(subdomain);
199 :
200 1333857 : for (const auto elem : _mesh.getMesh().active_local_subdomain_elements_ptr_range(subdomain))
201 : {
202 665348 : const auto elem_type = getElementTypeID(elem);
203 665348 : const auto eid = getContiguousElementID(elem);
204 :
205 665348 : _elem_info[eid].type = elem_type;
206 665348 : _elem_info[eid].id = eid;
207 665348 : _elem_info[eid].subdomain = sid;
208 :
209 2339212 : for (const auto i : make_range(_mesh.getMesh().n_elem_integers()))
210 1673864 : _extra_elem_ids(eid, i) = elem->get_extra_integer(i);
211 3161 : }
212 :
213 3161 : _starting_elem_id[sid] = *getSubdomainContiguousElementIDRange(subdomain).begin();
214 : }
215 :
216 : // Cache ElementInfo for one-layer ghost neighbors so side-based computations can identify
217 : // off-process neighbor element subdomains.
218 2492 : for (const auto & [ghost_elem, ghost_eid] : _maps->ghost_elem_id_mapping)
219 : {
220 32 : _elem_info[ghost_eid].type = getElementTypeID(ghost_elem);
221 32 : _elem_info[ghost_eid].id = ghost_eid;
222 32 : _elem_info[ghost_eid].subdomain = getContiguousSubdomainID(ghost_elem->subdomain_id());
223 : }
224 :
225 2460 : _elem_info.moveToDevice();
226 2460 : _extra_elem_ids.moveToDevice();
227 2460 : _starting_elem_id.moveToDevice();
228 :
229 : // Cache node data
230 :
231 2460 : const auto num_nodes = getNumLocalNodes();
232 :
233 2460 : _points.create(num_nodes);
234 2460 : _nodes.create(_mesh.getMaxNodesPerElem(), num_local_plus_one_neighbor_layer_elems);
235 2460 : _boundary_nodes.create(_mesh.meshBoundaryIds().size());
236 :
237 734666 : for (const auto node : getLocalNodes())
238 732206 : _points[getContiguousNodeID(node)] = *node;
239 :
240 1333156 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
241 : {
242 665348 : const auto eid = getContiguousElementID(elem);
243 :
244 3357034 : for (unsigned int node = 0; node < elem->n_nodes(); ++node)
245 2691686 : _nodes(node, eid) = getContiguousNodeID(elem->node_ptr(node));
246 2460 : }
247 :
248 12117 : for (const auto boundary : _mesh.meshBoundaryIds())
249 9657 : _boundary_nodes[getContiguousBoundaryID(boundary)].copySet<false, true>(
250 5257 : _maps->boundary_node_ids[boundary]);
251 :
252 2460 : _points.moveToDevice();
253 2460 : _nodes.moveToDevice();
254 2460 : _boundary_nodes.copyToDevice();
255 2460 : }
256 :
257 : ContiguousSubdomainID
258 1041199 : Mesh::getContiguousSubdomainID(const SubdomainID subdomain) const
259 : {
260 1041199 : return libmesh_map_find(_maps->subdomain_id_mapping, subdomain);
261 : }
262 :
263 : ContiguousBoundaryID
264 9897 : Mesh::getContiguousBoundaryID(const BoundaryID boundary) const
265 : {
266 9897 : return libmesh_map_find(_maps->boundary_id_mapping, boundary);
267 : }
268 :
269 : unsigned int
270 665380 : Mesh::getElementTypeID(const Elem * elem) const
271 : {
272 : mooseAssert(elem, "Element pointer is null");
273 :
274 665380 : return libmesh_map_find(_maps->elem_type_id_mapping, elem->type());
275 : }
276 :
277 : ContiguousElementID
278 4065436 : Mesh::getContiguousElementID(const Elem * elem) const
279 : {
280 : mooseAssert(elem, "Element pointer is null");
281 :
282 4065436 : if (elem->processor_id() == _mesh.processor_id())
283 4065404 : return elem->get_extra_integer(_elem_id_integer);
284 :
285 32 : return libmesh_map_find(_maps->ghost_elem_id_mapping, elem);
286 : }
287 :
288 : ContiguousNodeID
289 4999101 : Mesh::getContiguousNodeID(const Node * node) const
290 : {
291 : mooseAssert(node, "Node pointer is null");
292 :
293 4999101 : if (node->processor_id() == _mesh.processor_id())
294 4977731 : return node->get_extra_integer(_node_id_integer);
295 : else
296 21370 : return libmesh_map_find(_maps->ghost_node_id_mapping, node);
297 : }
298 :
299 : void
300 181 : Mesh::initElementGeometry()
301 : {
302 181 : const auto num_local_elems = getNumLocalElements();
303 :
304 181 : _elem_volume.create(num_local_elems);
305 181 : _elem_centroid.create(num_local_elems);
306 :
307 181 : _elem_volume = 0;
308 181 : _elem_centroid = Real3();
309 :
310 42927 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
311 : {
312 21373 : const auto eid = getContiguousElementID(elem);
313 21373 : const auto elem_centroid = elem->vertex_average();
314 : Real coord_factor;
315 21373 : ::coordTransformFactor(_mesh, elem->subdomain_id(), elem_centroid, coord_factor);
316 :
317 21373 : _elem_volume[eid] = elem->volume() * coord_factor;
318 21373 : _elem_centroid[eid] = elem_centroid;
319 181 : }
320 :
321 181 : _elem_volume.moveToDevice();
322 181 : _elem_centroid.moveToDevice();
323 :
324 181 : _element_geometry_initialized = true;
325 181 : }
326 :
327 : void
328 181 : Mesh::initElementSideGeometry()
329 : {
330 : mooseAssert(_element_geometry_initialized,
331 : "initElementGeometry() must be called before initElementSideGeometry().");
332 :
333 181 : const auto num_elems = getNumLocalElements();
334 181 : const auto max_sides = _mesh.getMaxSidesPerElem();
335 :
336 181 : _side_area.create(max_sides, num_elems);
337 181 : _side_centroid.create(max_sides, num_elems);
338 181 : _side_normal.create(max_sides, num_elems);
339 181 : _elem_centroid_to_side_centroid.create(max_sides, num_elems);
340 181 : _elem_centroid_to_side_centroid_distance.create(max_sides, num_elems);
341 181 : _elem_centroid_to_neighbor_centroid.create(max_sides, num_elems);
342 181 : _elem_centroid_to_neighbor_centroid_distance.create(max_sides, num_elems);
343 181 : _side_boundary_id.create(max_sides, num_elems);
344 :
345 181 : _side_area = 0;
346 181 : _side_centroid = Real3();
347 181 : _side_normal = Real3();
348 181 : _elem_centroid_to_side_centroid = Real3();
349 181 : _elem_centroid_to_side_centroid_distance = 0;
350 181 : _elem_centroid_to_neighbor_centroid = Real3();
351 181 : _elem_centroid_to_neighbor_centroid_distance = 0;
352 181 : _side_boundary_id = Moose::INVALID_BOUNDARY_ID;
353 :
354 181 : libMesh::ElemSideBuilder side_builder;
355 :
356 42927 : for (const auto elem : _mesh.getMesh().active_local_element_ptr_range())
357 : {
358 21373 : const auto eid = getContiguousElementID(elem);
359 21373 : const auto elem_centroid = elem->vertex_average();
360 :
361 102487 : for (unsigned int side = 0; side < elem->n_sides(); ++side)
362 : {
363 81114 : auto & face = side_builder(*elem, side);
364 81114 : const auto face_centroid = face.vertex_average();
365 81114 : const auto * const neighbor = elem->neighbor_ptr(side);
366 : Real coord_factor;
367 : mooseAssert(!neighbor || neighbor != libMesh::remote_elem,
368 : "First layer neighbors should be ghosted");
369 81114 : ::coordTransformFactor(_mesh,
370 40588 : elem->subdomain_id(),
371 : face_centroid,
372 : coord_factor,
373 38804 : neighbor ? neighbor->subdomain_id()
374 : : libMesh::Elem::invalid_subdomain_id);
375 :
376 81114 : ContiguousElementID neighbor_eid = libMesh::DofObject::invalid_id;
377 81114 : if (neighbor)
378 : {
379 77562 : neighbor_eid = getContiguousElementID(neighbor);
380 : mooseAssert(
381 : neighbor_eid != libMesh::DofObject::invalid_id,
382 : "We should have determined a contiguous element ID for our local element neighbors");
383 77562 : _elem_neighbor(side, eid) = neighbor_eid;
384 : }
385 :
386 81114 : const auto elem_centroid_to_side_centroid = face_centroid - elem_centroid;
387 40526 : const auto elem_centroid_to_neighbor_centroid =
388 40588 : neighbor ? neighbor->vertex_average() - elem_centroid : Real3();
389 :
390 81114 : _side_area(side, eid) = face.volume() * coord_factor;
391 81114 : _side_centroid(side, eid) = face_centroid;
392 81114 : _side_normal(side, eid) = elem->side_vertex_average_normal(side);
393 81114 : _elem_centroid_to_side_centroid(side, eid) = elem_centroid_to_side_centroid;
394 81114 : _elem_centroid_to_side_centroid_distance(side, eid) = elem_centroid_to_side_centroid.norm();
395 81114 : _elem_centroid_to_neighbor_centroid(side, eid) = elem_centroid_to_neighbor_centroid;
396 81114 : _elem_centroid_to_neighbor_centroid_distance(side, eid) =
397 40588 : neighbor ? elem_centroid_to_neighbor_centroid.norm() : Real(0);
398 :
399 81114 : const auto boundary_ids = _mesh.getBoundaryIDs(elem, side);
400 81114 : if (boundary_ids.size() > 1)
401 0 : mooseError("Kokkos element-side geometry does not support multiple boundary IDs on a "
402 : "single side. "
403 : "Element ",
404 : eid,
405 : " side ",
406 : side,
407 : " has ",
408 0 : boundary_ids.size(),
409 : " boundary IDs.");
410 81114 : if (!boundary_ids.empty())
411 4048 : _side_boundary_id(side, eid) = boundary_ids.front();
412 81114 : }
413 181 : }
414 :
415 181 : _elem_neighbor.moveToDevice();
416 181 : _side_area.moveToDevice();
417 181 : _side_centroid.moveToDevice();
418 181 : _side_normal.moveToDevice();
419 181 : _elem_centroid_to_side_centroid.moveToDevice();
420 181 : _elem_centroid_to_side_centroid_distance.moveToDevice();
421 181 : _elem_centroid_to_neighbor_centroid.moveToDevice();
422 181 : _elem_centroid_to_neighbor_centroid_distance.moveToDevice();
423 181 : _side_boundary_id.moveToDevice();
424 :
425 181 : _element_side_geometry_initialized = true;
426 181 : }
427 :
428 : } // namespace Moose::Kokkos
|
