DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
utils.h
Go to the documentation of this file.
1// Copyright (C) 2019-2024 Garth N. Wells
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include "Mesh.h"
10#include "Topology.h"
11#include "graphbuild.h"
12#include <algorithm>
13#include <basix/mdspan.hpp>
14#include <concepts>
15#include <dolfinx/graph/AdjacencyList.h>
16#include <dolfinx/graph/ordering.h>
17#include <dolfinx/graph/partition.h>
18#include <functional>
19#include <mpi.h>
20#include <numeric>
21#include <span>
22
25
26namespace dolfinx::fem
27{
29}
30
31namespace dolfinx::mesh
32{
33enum class CellType;
34
36enum class GhostMode : int
37{
38 none,
39 shared_facet,
40 shared_vertex
41};
42
43namespace impl
44{
50template <typename T>
51void reorder_list(std::span<T> list, std::span<const std::int32_t> nodemap)
52{
53 if (nodemap.empty())
54 return;
55
56 assert(list.size() % nodemap.size() == 0);
57 std::size_t degree = list.size() / nodemap.size();
58 const std::vector<T> orig(list.begin(), list.end());
59 for (std::size_t n = 0; n < nodemap.size(); ++n)
60 {
61 std::span links_old(orig.data() + n * degree, degree);
62 auto links_new = list.subspan(nodemap[n] * degree, degree);
63 std::ranges::copy(links_old, links_new.begin());
64 }
65}
66
80template <std::floating_point T>
81std::tuple<std::vector<std::int32_t>, std::vector<T>, std::vector<std::int32_t>>
83 std::span<const std::int32_t> facets)
84{
85 auto topology = mesh.topology();
86 assert(topology);
87 const int tdim = topology->dim();
88 if (dim == tdim)
89 {
90 throw std::runtime_error(
91 "Cannot use mesh::locate_entities_boundary (boundary) for cells.");
92 }
93
94 // Build set of vertices on boundary and set of boundary entities
95 mesh.topology_mutable()->create_connectivity(tdim - 1, 0);
96 mesh.topology_mutable()->create_connectivity(tdim - 1, dim);
97 std::vector<std::int32_t> vertices, entities;
98 {
99 auto f_to_v = topology->connectivity(tdim - 1, 0);
100 assert(f_to_v);
101 auto f_to_e = topology->connectivity(tdim - 1, dim);
102 assert(f_to_e);
103 for (auto f : facets)
104 {
105 auto v = f_to_v->links(f);
106 vertices.insert(vertices.end(), v.begin(), v.end());
107 auto e = f_to_e->links(f);
108 entities.insert(entities.end(), e.begin(), e.end());
109 }
110
111 // Build vector of boundary vertices
112 {
113 std::ranges::sort(vertices);
114 auto [unique_end, range_end] = std::ranges::unique(vertices);
115 vertices.erase(unique_end, range_end);
116 }
117
118 {
119 std::ranges::sort(entities);
120 auto [unique_end, range_end] = std::ranges::unique(entities);
121 entities.erase(unique_end, range_end);
122 }
123 }
124
125 // Get geometry data
126 auto x_dofmap = mesh.geometry().dofmap();
127 std::span<const T> x_nodes = mesh.geometry().x();
128
129 // Get all vertex 'node' indices
130 mesh.topology_mutable()->create_connectivity(0, tdim);
131 mesh.topology_mutable()->create_connectivity(tdim, 0);
132 auto v_to_c = topology->connectivity(0, tdim);
133 assert(v_to_c);
134 auto c_to_v = topology->connectivity(tdim, 0);
135 assert(c_to_v);
136 std::vector<T> x_vertices(3 * vertices.size(), -1.0);
137 std::vector<std::int32_t> vertex_to_pos(v_to_c->num_nodes(), -1);
138 for (std::size_t i = 0; i < vertices.size(); ++i)
139 {
140 const std::int32_t v = vertices[i];
141
142 // Get first cell and find position
143 const std::int32_t c = v_to_c->links(v).front();
144 auto cell_vertices = c_to_v->links(c);
145 auto it = std::find(cell_vertices.begin(), cell_vertices.end(), v);
146 assert(it != cell_vertices.end());
147 const std::size_t local_pos = std::distance(cell_vertices.begin(), it);
148
149 auto dofs = md::submdspan(x_dofmap, c, md::full_extent);
150 for (std::size_t j = 0; j < 3; ++j)
151 x_vertices[j * vertices.size() + i] = x_nodes[3 * dofs[local_pos] + j];
152 vertex_to_pos[v] = i;
153 }
154
155 return {std::move(entities), std::move(x_vertices), std::move(vertex_to_pos)};
156}
157
158} // namespace impl
159
173std::vector<std::int32_t> exterior_facet_indices(const Topology& topology,
174 int facet_type_idx);
175
187std::vector<std::int32_t> exterior_facet_indices(const Topology& topology);
188
206using CellPartitionFunction = std::function<graph::AdjacencyList<std::int32_t>(
207 MPI_Comm comm, int nparts, const std::vector<CellType>& cell_types,
208 const std::vector<std::span<const std::int64_t>>& cells)>;
209
214using CellReorderFunction = std::function<std::vector<std::int32_t>(
216
226std::vector<std::int64_t> extract_topology(CellType cell_type,
227 const fem::ElementDofLayout& layout,
228 std::span<const std::int64_t> cells);
229
238template <std::floating_point T>
239std::vector<T> h(const Mesh<T>& mesh, std::span<const std::int32_t> entities,
240 int dim)
241{
242 if (entities.empty())
243 return std::vector<T>();
244 if (dim == 0)
245 return std::vector<T>(entities.size(), 0);
246
247 // Get the geometry dofs for the vertices of each entity
248 const auto [vertex_xdofs, xdof_shape]
249 = entities_to_geometry(mesh, dim, entities, false);
250
251 // Get the geometry coordinate
252 std::span<const T> x = mesh.geometry().x();
253
254 // Function to compute the length of (p0 - p1)
255 auto delta_norm = [](auto&& p0, auto&& p1)
256 {
257 T norm = 0;
258 for (std::size_t i = 0; i < 3; ++i)
259 norm += (p0[i] - p1[i]) * (p0[i] - p1[i]);
260 return std::sqrt(norm);
261 };
262
263 // Compute greatest distance between any to vertices
264 assert(dim > 0);
265 std::vector<T> h(entities.size(), 0);
266 for (std::size_t e = 0; e < entities.size(); ++e)
267 {
268 // Get geometry 'dof' for each vertex of entity e
269 std::span<const std::int32_t> e_vertices(
270 vertex_xdofs.data() + e * xdof_shape[1], xdof_shape[1]);
271
272 // Compute maximum distance between any two vertices
273 for (std::size_t i = 0; i < e_vertices.size(); ++i)
274 {
275 std::span<const T, 3> p0(x.data() + 3 * e_vertices[i], 3);
276 for (std::size_t j = i + 1; j < e_vertices.size(); ++j)
277 {
278 std::span<const T, 3> p1(x.data() + 3 * e_vertices[j], 3);
279 h[e] = std::max(h[e], delta_norm(p0, p1));
280 }
281 }
282 }
283
284 return h;
285}
286
290template <std::floating_point T>
291std::vector<T> cell_normals(const Mesh<T>& mesh, int dim,
292 std::span<const std::int32_t> entities)
293{
294 if (entities.empty())
295 return std::vector<T>();
296
297 auto topology = mesh.topology();
298 assert(topology);
299 if (topology->cell_type() == CellType::prism and dim == 2)
300 {
301 throw std::runtime_error(
302 "Cell normal computation for prism cells not yet supported.");
303 }
304
305 const int gdim = mesh.geometry().dim();
306 const CellType type = cell_entity_type(topology->cell_type(), dim, 0);
307
308 // Find geometry nodes for topology entities
309 std::span<const T> x = mesh.geometry().x();
310 const auto [geometry_entities, eshape]
311 = entities_to_geometry(mesh, dim, entities, false);
312
313 std::vector<T> n(entities.size() * 3);
314 switch (type)
315 {
316 case CellType::interval:
317 {
318 if (gdim > 2)
319 throw std::invalid_argument("Interval cell normal undefined in 3D.");
320 for (std::size_t i = 0; i < entities.size(); ++i)
321 {
322 // Get the two vertices as points
323 std::array vertices{geometry_entities[i * eshape[1]],
324 geometry_entities[i * eshape[1] + 1]};
325 std::array p = {std::span<const T, 3>(x.data() + 3 * vertices[0], 3),
326 std::span<const T, 3>(x.data() + 3 * vertices[1], 3)};
327
328 // Define normal by rotating tangent counter-clockwise
329 std::array<T, 3> t;
330 std::ranges::transform(p[1], p[0], t.begin(),
331 [](auto x, auto y) { return x - y; });
332
333 T norm = std::sqrt(t[0] * t[0] + t[1] * t[1]);
334 std::span<T, 3> ni(n.data() + 3 * i, 3);
335 ni[0] = -t[1] / norm;
336 ni[1] = t[0] / norm;
337 ni[2] = 0.0;
338 }
339 return n;
340 }
341 case CellType::triangle:
342 {
343 for (std::size_t i = 0; i < entities.size(); ++i)
344 {
345 // Get the three vertices as points
346 std::array vertices = {geometry_entities[i * eshape[1] + 0],
347 geometry_entities[i * eshape[1] + 1],
348 geometry_entities[i * eshape[1] + 2]};
349 std::array p = {std::span<const T, 3>(x.data() + 3 * vertices[0], 3),
350 std::span<const T, 3>(x.data() + 3 * vertices[1], 3),
351 std::span<const T, 3>(x.data() + 3 * vertices[2], 3)};
352
353 // Compute (p1 - p0) and (p2 - p0)
354 std::array<T, 3> dp1, dp2;
355 std::ranges::transform(p[1], p[0], dp1.begin(),
356 [](auto x, auto y) { return x - y; });
357 std::ranges::transform(p[2], p[0], dp2.begin(),
358 [](auto x, auto y) { return x - y; });
359
360 // Define cell normal via cross product of first two edges
361 std::array<T, 3> ni = math::cross(dp1, dp2);
362 T norm = std::sqrt(ni[0] * ni[0] + ni[1] * ni[1] + ni[2] * ni[2]);
363 std::ranges::transform(ni, std::next(n.begin(), 3 * i),
364 [norm](auto x) { return x / norm; });
365 }
366
367 return n;
368 }
369 case CellType::quadrilateral:
370 {
371 // TODO: check
372 for (std::size_t i = 0; i < entities.size(); ++i)
373 {
374 // Get the three vertices as points
375 std::array vertices = {geometry_entities[i * eshape[1] + 0],
376 geometry_entities[i * eshape[1] + 1],
377 geometry_entities[i * eshape[1] + 2]};
378 std::array p = {std::span<const T, 3>(x.data() + 3 * vertices[0], 3),
379 std::span<const T, 3>(x.data() + 3 * vertices[1], 3),
380 std::span<const T, 3>(x.data() + 3 * vertices[2], 3)};
381
382 // Compute (p1 - p0) and (p2 - p0)
383 std::array<T, 3> dp1, dp2;
384 std::ranges::transform(p[1], p[0], dp1.begin(),
385 [](auto x, auto y) { return x - y; });
386 std::ranges::transform(p[2], p[0], dp2.begin(),
387 [](auto x, auto y) { return x - y; });
388
389 // Define cell normal via cross product of first two edges
390 std::array<T, 3> ni = math::cross(dp1, dp2);
391 T norm = std::sqrt(ni[0] * ni[0] + ni[1] * ni[1] + ni[2] * ni[2]);
392 std::ranges::transform(ni, std::next(n.begin(), 3 * i),
393 [norm](auto x) { return x / norm; });
394 }
395
396 return n;
397 }
398 default:
399 throw std::invalid_argument(
400 "cell_normal not supported for this cell type.");
401 }
402}
403
407template <std::floating_point T>
408std::vector<T> compute_midpoints(const Mesh<T>& mesh, int dim,
409 std::span<const std::int32_t> entities)
410{
411 if (entities.empty())
412 return std::vector<T>();
413
414 std::span<const T> x = mesh.geometry().x();
415
416 // Build map from entity -> geometry dof
417 const auto [e_to_g, eshape]
418 = entities_to_geometry(mesh, dim, entities, false);
419
420 std::vector<T> x_mid(entities.size() * 3, 0);
421 for (std::size_t e = 0; e < entities.size(); ++e)
422 {
423 std::span<T, 3> p(x_mid.data() + 3 * e, 3);
424 std::span<const std::int32_t> rows(e_to_g.data() + e * eshape[1],
425 eshape[1]);
426 for (auto row : rows)
427 {
428 std::span<const T, 3> xg(x.data() + 3 * row, 3);
429 std::ranges::transform(p, xg, p.begin(),
430 [size = rows.size()](auto x, auto y)
431 { return x + y / size; });
432 }
433 }
434
435 return x_mid;
436}
437
438namespace impl
439{
444template <std::floating_point T>
445std::pair<std::vector<T>, std::array<std::size_t, 2>>
447{
448 auto topology = mesh.topology();
449 assert(topology);
450 const int tdim = topology->dim();
451
452 // Create entities and connectivities
453
454 // Get all vertex 'node' indices
455 const std::int32_t num_vertices = topology->index_map(0)->size_local()
456 + topology->index_map(0)->num_ghosts();
457
458 std::vector<std::int32_t> vertex_to_node(num_vertices);
459 for (int cell_type_idx = 0,
460 num_cell_types = topology->entity_types(tdim).size();
461 cell_type_idx < num_cell_types; ++cell_type_idx)
462 {
463 auto x_dofmap = mesh.geometry().dofmap(cell_type_idx);
464 auto c_to_v = topology->connectivity({tdim, cell_type_idx}, {0, 0});
465 assert(c_to_v);
466 for (int c = 0; c < c_to_v->num_nodes(); ++c)
467 {
468 auto x_dofs = md::submdspan(x_dofmap, c, md::full_extent);
469 auto vertices = c_to_v->links(c);
470 for (std::size_t i = 0; i < vertices.size(); ++i)
471 vertex_to_node[vertices[i]] = x_dofs[i];
472 }
473 }
474
475 // Pack coordinates of vertices
476 std::span<const T> x_nodes = mesh.geometry().x();
477 std::vector<T> x_vertices(3 * vertex_to_node.size(), 0.0);
478 for (std::size_t i = 0; i < vertex_to_node.size(); ++i)
479 {
480 std::int32_t pos = 3 * vertex_to_node[i];
481 for (std::size_t j = 0; j < 3; ++j)
482 x_vertices[j * vertex_to_node.size() + i] = x_nodes[pos + j];
483 }
484
485 return {std::move(x_vertices), {3, vertex_to_node.size()}};
486}
487
488} // namespace impl
489
491template <typename Fn, typename T>
492concept MarkerFn = std::is_invocable_r<
493 std::vector<std::int8_t>, Fn,
494 md::mdspan<const T,
495 md::extents<std::size_t, 3, md::dynamic_extent>>>::value;
496
512template <std::floating_point T, MarkerFn<T> U>
513std::vector<std::int32_t> locate_entities(const Mesh<T>& mesh, int dim,
514 U marker, int entity_type_idx)
515{
516
517 using cmdspan3x_t
518 = md::mdspan<const T, md::extents<std::size_t, 3, md::dynamic_extent>>;
519
520 // Run marker function on vertex coordinates
521 const auto [xdata, xshape] = impl::compute_vertex_coords(mesh);
522
523 cmdspan3x_t x(xdata.data(), xshape);
524 const std::vector<std::int8_t> marked = marker(x);
525 if (marked.size() != x.extent(1))
526 throw std::runtime_error("Length of array of markers is wrong.");
527
528 auto topology = mesh.topology();
529 assert(topology);
530 const int tdim = topology->dim();
531
532 mesh.topology_mutable()->create_entities(dim);
533 if (dim < tdim)
534 mesh.topology_mutable()->create_connectivity(dim, 0);
535
536 // Iterate over entities of dimension 'dim' to build vector of marked
537 // entities
538 auto e_to_v = topology->connectivity({dim, entity_type_idx}, {0, 0});
539 assert(e_to_v);
540 std::vector<std::int32_t> entities;
541 for (int e = 0; e < e_to_v->num_nodes(); ++e)
542 {
543 // Iterate over entity vertices
544 bool all_vertices_marked = true;
545 for (std::int32_t v : e_to_v->links(e))
546 {
547 if (!marked[v])
548 {
549 all_vertices_marked = false;
550 break;
551 }
552 }
553
554 if (all_vertices_marked)
555 entities.push_back(e);
556 }
557
558 return entities;
559}
560
574template <std::floating_point T, MarkerFn<T> U>
575std::vector<std::int32_t> locate_entities(const Mesh<T>& mesh, int dim,
576 U marker)
577{
578 const int num_entity_types = mesh.topology()->entity_types(dim).size();
579 if (num_entity_types > 1)
580 {
581 throw std::runtime_error(
582 "Multiple entity types of this dimension. Specify entity type index");
583 }
584 return locate_entities(mesh, dim, marker, 0);
585}
586
610template <std::floating_point T, MarkerFn<T> U>
611std::vector<std::int32_t> locate_entities_boundary(const Mesh<T>& mesh, int dim,
612 U marker)
613{
614 // TODO Rewrite this function, it should be possible to simplify considerably
615 auto topology = mesh.topology();
616 assert(topology);
617 int tdim = topology->dim();
618 if (dim == tdim)
619 {
620 throw std::runtime_error(
621 "Cannot use mesh::locate_entities_boundary (boundary) for cells.");
622 }
623
624 // Compute list of boundary facets
625 mesh.topology_mutable()->create_entities(tdim - 1);
626 mesh.topology_mutable()->create_connectivity(tdim - 1, tdim);
627 std::vector<std::int32_t> boundary_facets = exterior_facet_indices(*topology);
628
629 using cmdspan3x_t
630 = md::mdspan<const T, md::extents<std::size_t, 3, md::dynamic_extent>>;
631
632 // Run marker function on the vertex coordinates
633 auto [facet_entities, xdata, vertex_to_pos]
634 = impl::compute_vertex_coords_boundary(mesh, dim, boundary_facets);
635 cmdspan3x_t x(xdata.data(), 3, xdata.size() / 3);
636 std::vector<std::int8_t> marked = marker(x);
637 if (marked.size() != x.extent(1))
638 throw std::runtime_error("Length of array of markers is wrong.");
639
640 // Loop over entities and check vertex markers
641 mesh.topology_mutable()->create_entities(dim);
642 auto e_to_v = topology->connectivity(dim, 0);
643 assert(e_to_v);
644 std::vector<std::int32_t> entities;
645 for (auto e : facet_entities)
646 {
647 // Iterate over entity vertices
648 bool all_vertices_marked = true;
649 for (auto v : e_to_v->links(e))
650 {
651 const std::int32_t pos = vertex_to_pos[v];
652 if (!marked[pos])
653 {
654 all_vertices_marked = false;
655 break;
656 }
657 }
658
659 // Mark facet with all vertices marked
660 if (all_vertices_marked)
661 entities.push_back(e);
662 }
663
664 return entities;
665}
666
685template <std::floating_point T>
686std::pair<std::vector<std::int32_t>, std::array<std::size_t, 2>>
688 std::span<const std::int32_t> entities,
689 bool permute = false)
690{
691 auto topology = mesh.topology();
692 assert(topology);
693 CellType cell_type = topology->cell_type();
694 if (cell_type == CellType::prism and dim == 2)
695 {
696 throw std::runtime_error(
697 "mesh::entities_to_geometry for prism cells not yet supported.");
698 }
699
700 const int tdim = topology->dim();
701 const Geometry<T>& geometry = mesh.geometry();
702 auto xdofs = geometry.dofmap();
703
704 // Get the DOF layout and the number of DOFs per entity
705 const fem::CoordinateElement<T>& coord_ele = geometry.cmap();
706 const fem::ElementDofLayout layout = coord_ele.create_dof_layout();
707 const std::size_t num_entity_dofs = layout.num_entity_closure_dofs(dim);
708 std::vector<std::int32_t> entity_xdofs;
709 entity_xdofs.reserve(entities.size() * num_entity_dofs);
710 std::array<std::size_t, 2> eshape{entities.size(), num_entity_dofs};
711
712 // Get the element's closure DOFs
713 const std::vector<std::vector<std::vector<int>>>& closure_dofs_all
714 = layout.entity_closure_dofs_all();
715
716 // Special case when dim == tdim (cells)
717 if (dim == tdim)
718 {
719 for (std::int32_t c : entities)
720 {
721 // Extract degrees of freedom
722 auto x_c = md::submdspan(xdofs, c, md::full_extent);
723 for (std::int32_t entity_dof : closure_dofs_all[tdim][0])
724 entity_xdofs.push_back(x_c[entity_dof]);
725 }
726
727 return {std::move(entity_xdofs), eshape};
728 }
729
730 assert(dim != tdim);
731
732 auto e_to_c = topology->connectivity(dim, tdim);
733 if (!e_to_c)
734 {
735 throw std::runtime_error(
736 "Entity-to-cell connectivity has not been computed. Missing dims "
737 + std::to_string(dim) + "->" + std::to_string(tdim));
738 }
739
740 auto c_to_e = topology->connectivity(tdim, dim);
741 if (!c_to_e)
742 {
743 throw std::runtime_error(
744 "Cell-to-entity connectivity has not been computed. Missing dims "
745 + std::to_string(tdim) + "->" + std::to_string(dim));
746 }
747
748 // Get the cell info, which is needed to permute the closure dofs
749 std::span<const std::uint32_t> cell_info;
750 if (permute)
751 cell_info = std::span(mesh.topology()->get_cell_permutation_info());
752
753 for (std::int32_t e : entities)
754 {
755 // Get a cell connected to the entity
756 assert(!e_to_c->links(e).empty());
757 std::int32_t c = e_to_c->links(e).front();
758
759 // Get the local index of the entity
760 std::span<const std::int32_t> cell_entities = c_to_e->links(c);
761 auto it = std::find(cell_entities.begin(), cell_entities.end(), e);
762 assert(it != cell_entities.end());
763 std::size_t local_entity = std::distance(cell_entities.begin(), it);
764
765 // Cell sub-entities must be permuted so that their local
766 // orientation agrees with their global orientation
767 std::vector<std::int32_t> closure_dofs(closure_dofs_all[dim][local_entity]);
768 if (permute)
769 {
770 mesh::CellType entity_type
771 = mesh::cell_entity_type(cell_type, dim, local_entity);
772 coord_ele.permute_subentity_closure(closure_dofs, cell_info[c],
773 entity_type, local_entity);
774 }
775
776 // Extract degrees of freedom
777 auto x_c = md::submdspan(xdofs, c, md::full_extent);
778 for (std::int32_t entity_dof : closure_dofs)
779 entity_xdofs.push_back(x_c[entity_dof]);
780 }
781
782 return {std::move(entity_xdofs), eshape};
783}
784
790 = mesh::GhostMode::none,
791 const graph::partition_fn& partfn
793
801std::vector<std::int32_t>
802compute_incident_entities(const Topology& topology,
803 std::span<const std::int32_t> entities, int d0,
804 int d1);
805
847template <typename U>
849 MPI_Comm comm, MPI_Comm commt,
850 std::vector<std::span<const std::int64_t>> cells,
851 const std::vector<fem::CoordinateElement<
852 typename std::remove_reference_t<typename U::value_type>>>& elements,
853 MPI_Comm commg, const U& x, std::array<std::size_t, 2> xshape,
854 const CellPartitionFunction& partitioner,
855 const CellReorderFunction& reorder_fn = graph::reorder_gps)
856{
857 assert(cells.size() == elements.size());
858 std::vector<CellType> celltypes;
859 std::ranges::transform(elements, std::back_inserter(celltypes),
860 [](auto e) { return e.cell_shape(); });
861 std::vector<fem::ElementDofLayout> doflayouts;
862 std::ranges::transform(elements, std::back_inserter(doflayouts),
863 [](auto e) { return e.create_dof_layout(); });
864
865 // Note: `extract_topology` extracts topology data, i.e. just the
866 // vertices. For P1 geometry this should just be the identity
867 // operator. For other elements the filtered lists may have 'gaps',
868 // i.e. the indices might not be contiguous.
869 //
870 // `extract_topology` could be skipped for 'P1 geometry' elements
871
872 std::int32_t num_cell_types = cells.size();
873
874 // -- Partition topology across ranks of comm
875 std::vector<std::vector<std::int64_t>> cells1(num_cell_types);
876 std::vector<std::vector<std::int64_t>> original_idx1(num_cell_types);
877 std::vector<std::vector<int>> ghost_owners(num_cell_types);
878 if (partitioner)
879 {
880 spdlog::info("Using partitioner with cell data ({} cell types)",
881 num_cell_types);
883 if (commt != MPI_COMM_NULL)
884 {
885 int size = dolfinx::MPI::size(comm);
886 std::vector<std::vector<std::int64_t>> t(num_cell_types);
887 std::vector<std::span<const std::int64_t>> tspan(num_cell_types);
888 for (std::int32_t i = 0; i < num_cell_types; ++i)
889 {
890 t[i] = extract_topology(celltypes[i], doflayouts[i], cells[i]);
891 tspan[i] = std::span(t[i]);
892 }
893 dest = partitioner(commt, size, celltypes, tspan);
894 }
895
896 std::int32_t cell_offset = 0;
897 for (std::int32_t i = 0; i < num_cell_types; ++i)
898 {
899 std::size_t num_cell_nodes = doflayouts[i].num_dofs();
900 assert(cells[i].size() % num_cell_nodes == 0);
901 std::size_t num_cells = cells[i].size() / num_cell_nodes;
902
903 // Extract destination AdjacencyList for this cell type
904 std::vector<std::int32_t> offsets_i(
905 std::next(dest.offsets().begin(), cell_offset),
906 std::next(dest.offsets().begin(), cell_offset + num_cells + 1));
907 std::vector<std::int32_t> data_i(
908 std::next(dest.array().begin(), offsets_i.front()),
909 std::next(dest.array().begin(), offsets_i.back()));
910 std::int32_t offset_0 = offsets_i.front();
911 std::ranges::for_each(offsets_i,
912 [&offset_0](std::int32_t& j) { j -= offset_0; });
913 graph::AdjacencyList<std::int32_t> dest_i(data_i, offsets_i);
914 cell_offset += num_cells;
915
916 // Distribute cells (topology, includes higher-order 'nodes') to
917 // destination rank
918 std::vector<int> src_ranks;
919 std::tie(cells1[i], src_ranks, original_idx1[i], ghost_owners[i])
920 = graph::build::distribute(comm, cells[i],
921 {num_cells, num_cell_nodes}, dest_i);
922 spdlog::debug("Got {} cells from distribution", cells1[i].size());
923 }
924 }
925 else
926 {
927 // No partitioning, construct a global index
928 std::int64_t num_owned = 0;
929 for (std::int32_t i = 0; i < num_cell_types; ++i)
930 {
931 cells1[i] = std::vector<std::int64_t>(cells[i].begin(), cells[i].end());
932 std::int32_t num_cell_nodes = doflayouts[i].num_dofs();
933 assert(cells1[i].size() % num_cell_nodes == 0);
934 original_idx1[i].resize(cells1[i].size() / num_cell_nodes);
935 num_owned += original_idx1[i].size();
936 }
937
938 // Add on global offset
939 std::int64_t global_offset = 0;
940 MPI_Exscan(&num_owned, &global_offset, 1, MPI_INT64_T, MPI_SUM, comm);
941 for (std::int32_t i = 0; i < num_cell_types; ++i)
942 {
943 std::iota(original_idx1[i].begin(), original_idx1[i].end(),
944 global_offset);
945 global_offset += original_idx1[i].size();
946 }
947 }
948
949 // Extract cell 'topology', i.e. extract the vertices for each cell
950 // and discard any 'higher-order' nodes
951 std::vector<std::vector<std::int64_t>> cells1_v(num_cell_types);
952 for (std::int32_t i = 0; i < num_cell_types; ++i)
953 {
954 cells1_v[i] = extract_topology(celltypes[i], doflayouts[i], cells1[i]);
955 spdlog::info("Extract basic topology: {}->{}", cells1[i].size(),
956 cells1_v[i].size());
957 }
958
959 // Build local dual graph for owned cells to (i) get list of vertices
960 // on the process boundary and (ii) apply re-ordering to cells for
961 // locality
962 auto boundary_v_fn = [](const std::vector<CellType>& celltypes,
963 const std::vector<fem::ElementDofLayout>& doflayouts,
964 const std::vector<std::vector<int>>& ghost_owners,
965 std::vector<std::vector<std::int64_t>>& cells1,
966 std::vector<std::vector<std::int64_t>>& cells1_v,
967 std::vector<std::vector<std::int64_t>>& original_idx1,
968 const CellReorderFunction& reorder_fn)
969 {
970 spdlog::info("Build local dual graphs, re-order cells, and compute process "
971 "boundary vertices.");
972
973 std::vector<std::pair<std::vector<std::int64_t>, int>> facets;
974
975 // Build lists of cells (by cell type) that excludes ghosts
976 std::vector<std::span<const std::int64_t>> cells1_v_local;
977 for (std::size_t i = 0; i < celltypes.size(); ++i)
978 {
979 int num_cell_vertices = mesh::num_cell_vertices(celltypes[i]);
980 std::size_t num_owned_cells
981 = cells1_v[i].size() / num_cell_vertices - ghost_owners[i].size();
982 cells1_v_local.emplace_back(cells1_v[i].data(),
983 num_owned_cells * num_cell_vertices);
984
985 // Build local dual graph for cell type
986 auto [graph, unmatched_facets, max_v, _facet_attached_cells]
987 = build_local_dual_graph(std::vector{celltypes[i]},
988 std::vector{cells1_v_local.back()});
989
990 // Store unmatched_facets for current cell type
991 facets.emplace_back(std::move(unmatched_facets), max_v);
992
993 // Compute re-ordering of graph
994 const std::vector<std::int32_t> remap = reorder_fn(graph);
995
996 // Update 'original' indices
997 const std::vector<std::int64_t>& orig_idx = original_idx1[i];
998 std::vector<std::int64_t> _original_idx(orig_idx.size());
999 std::copy_n(orig_idx.rbegin(), ghost_owners[i].size(),
1000 _original_idx.rbegin());
1001 {
1002 for (std::size_t j = 0; j < remap.size(); ++j)
1003 _original_idx[remap[j]] = orig_idx[j];
1004 }
1005 original_idx1[i] = _original_idx;
1006
1007 // Reorder cells
1009 std::span(cells1_v[i].data(), remap.size() * num_cell_vertices),
1010 remap);
1012 std::span(cells1[i].data(), remap.size() * doflayouts[i].num_dofs()),
1013 remap);
1014 }
1015
1016 if (facets.size() == 1) // Optimisation for single cell type
1017 {
1018 std::vector<std::int64_t>& vertices = facets.front().first;
1019
1020 // Remove duplicated vertex indices
1021 std::ranges::sort(vertices);
1022 auto [unique_end, range_end] = std::ranges::unique(vertices);
1023 vertices.erase(unique_end, range_end);
1024
1025 // Remove -1 if it appears as first entity. This can happen in
1026 // mixed topology meshes where '-1' is used to pad facet data when
1027 // cells facets have differing numbers of vertices.
1028 if (!vertices.empty() and vertices.front() == -1)
1029 vertices.erase(vertices.begin());
1030
1031 return vertices;
1032 }
1033 else
1034 {
1035 // Pack 'unmatched' facets for all cell types into single array
1036 // (facets0)
1037 std::vector<std::int64_t> facets0;
1038 facets0.reserve(std::accumulate(facets.begin(), facets.end(),
1039 std::size_t(0), [](std::size_t x, auto& y)
1040 { return x + y.first.size(); }));
1041 int max_v = std::ranges::max_element(facets, [](auto& a, auto& b)
1042 { return a.second < b.second; })
1043 ->second;
1044 for (const auto& [v_data, num_v] : facets)
1045 {
1046 for (auto it = v_data.begin(); it != v_data.end(); it += num_v)
1047 {
1048 facets0.insert(facets0.end(), it, std::next(it, num_v));
1049 facets0.insert(facets0.end(), max_v - num_v, -1);
1050 }
1051 }
1052
1053 // Compute row permutation
1054 const std::vector<std::int32_t> perm = dolfinx::sort_by_perm(
1055 std::span<const std::int64_t>(facets0), max_v);
1056
1057 // For facets in facets0 that appear only once, store the facet
1058 // vertices
1059 std::vector<std::int64_t> vertices;
1060 auto it = perm.begin();
1061 while (it != perm.end())
1062 {
1063 // Find iterator to next facet different from f and trim any -1
1064 // padding
1065 std::span _f(facets0.data() + (*it) * max_v, max_v);
1066 auto end = std::find_if(_f.rbegin(), _f.rend(),
1067 [](auto a) { return a >= 0; });
1068 auto f = _f.first(std::distance(end, _f.rend()));
1069
1070 auto it1 = std::find_if_not(
1071 it, perm.end(),
1072 [f, max_v, it0 = facets0.begin()](auto p) -> bool
1073 {
1074 return std::equal(f.begin(), f.end(), std::next(it0, p * max_v));
1075 });
1076
1077 // If no repeated facet found, insert f vertices
1078 if (std::distance(it, it1) == 1)
1079 vertices.insert(vertices.end(), f.begin(), f.end());
1080 else if (std::distance(it, it1) > 2)
1081 throw std::runtime_error("More than two matching facets found.");
1082
1083 // Advance iterator
1084 it = it1;
1085 }
1086
1087 // Remove duplicate indices
1088 std::ranges::sort(vertices);
1089 auto [unique_end, range_end] = std::ranges::unique(vertices);
1090 vertices.erase(unique_end, range_end);
1091
1092 return vertices;
1093 }
1094 };
1095
1096 const std::vector<std::int64_t> boundary_v
1097 = boundary_v_fn(celltypes, doflayouts, ghost_owners, cells1, cells1_v,
1098 original_idx1, reorder_fn);
1099
1100 spdlog::debug("Got {} boundary vertices", boundary_v.size());
1101
1102 // Create Topology
1103 std::vector<std::span<const std::int64_t>> cells1_v_span;
1104 std::ranges::transform(cells1_v, std::back_inserter(cells1_v_span),
1105 [](auto& c) { return std::span(c); });
1106 std::vector<std::span<const std::int64_t>> original_idx1_span;
1107 std::ranges::transform(original_idx1, std::back_inserter(original_idx1_span),
1108 [](auto& c) { return std::span(c); });
1109 std::vector<std::span<const int>> ghost_owners_span;
1110 std::ranges::transform(ghost_owners, std::back_inserter(ghost_owners_span),
1111 [](auto& c) { return std::span(c); });
1112 Topology topology
1113 = create_topology(comm, celltypes, cells1_v_span, original_idx1_span,
1114 ghost_owners_span, boundary_v);
1115
1116 // Create connectivities required higher-order geometries for creating
1117 // a Geometry object
1118 for (int i = 0; i < num_cell_types; ++i)
1119 {
1120 for (int e = 1; e < topology.dim(); ++e)
1121 {
1122 if (doflayouts[i].num_entity_dofs(e) > 0)
1123 topology.create_entities(e);
1124 }
1125
1126 if (elements[i].needs_dof_permutations())
1127 topology.create_entity_permutations();
1128 }
1129
1130 // Build list of unique (global) node indices from cells1 and
1131 // distribute coordinate data
1132 std::vector<std::int64_t> nodes1, nodes2;
1133 for (std::vector<std::int64_t>& c : cells1)
1134 nodes1.insert(nodes1.end(), c.begin(), c.end());
1135 for (std::vector<std::int64_t>& c : cells1)
1136 nodes2.insert(nodes2.end(), c.begin(), c.end());
1137
1138 dolfinx::radix_sort(nodes1);
1139 auto [unique_end, range_end] = std::ranges::unique(nodes1);
1140 nodes1.erase(unique_end, range_end);
1141
1142 std::vector coords
1143 = dolfinx::MPI::distribute_data(comm, nodes1, commg, x, xshape[1]);
1144
1145 // Create geometry object
1147 = create_geometry(topology, elements, nodes1, nodes2, coords, xshape[1]);
1148
1149 return Mesh(comm, std::make_shared<Topology>(std::move(topology)),
1150 std::move(geometry));
1151}
1152
1188template <typename U>
1190 MPI_Comm comm, MPI_Comm commt, std::span<const std::int64_t> cells,
1192 typename std::remove_reference_t<typename U::value_type>>& element,
1193 MPI_Comm commg, const U& x, std::array<std::size_t, 2> xshape,
1194 const CellPartitionFunction& partitioner,
1195 const CellReorderFunction& reorder_fn = graph::reorder_gps)
1196{
1197 return create_mesh(comm, commt, std::vector{cells}, std::vector{element},
1198 commg, x, xshape, partitioner, reorder_fn);
1199}
1200
1219template <typename U>
1220Mesh<typename std::remove_reference_t<typename U::value_type>>
1221create_mesh(MPI_Comm comm, std::span<const std::int64_t> cells,
1223 std::remove_reference_t<typename U::value_type>>& elements,
1224 const U& x, std::array<std::size_t, 2> xshape, GhostMode ghost_mode)
1225{
1226 if (dolfinx::MPI::size(comm) == 1)
1227 return create_mesh(comm, comm, std::vector{cells}, std::vector{elements},
1228 comm, x, xshape, nullptr);
1229 else
1230 {
1231 return create_mesh(comm, comm, std::vector{cells}, std::vector{elements},
1232 comm, x, xshape, create_cell_partitioner(ghost_mode));
1233 }
1234}
1235
1249template <std::floating_point T>
1250std::pair<Geometry<T>, std::vector<int32_t>>
1252 std::span<const std::int32_t> subentity_to_entity)
1253{
1254 const Geometry<T>& geometry = mesh.geometry();
1255
1256 // Get the geometry dofs in the sub-geometry based on the entities in
1257 // sub-geometry
1258 const fem::ElementDofLayout layout = geometry.cmap().create_dof_layout();
1259
1260 const std::vector<std::int32_t> x_indices
1261 = entities_to_geometry(mesh, dim, subentity_to_entity, true).first;
1262
1263 std::vector<std::int32_t> sub_x_dofs = x_indices;
1264 std::ranges::sort(sub_x_dofs);
1265 auto [unique_end, range_end] = std::ranges::unique(sub_x_dofs);
1266 sub_x_dofs.erase(unique_end, range_end);
1267
1268 // Get the sub-geometry dofs owned by this process
1269 auto x_index_map = geometry.index_map();
1270 assert(x_index_map);
1271
1272 std::shared_ptr<common::IndexMap> sub_x_dof_index_map;
1273 std::vector<std::int32_t> subx_to_x_dofmap;
1274 {
1275 auto [map, new_to_old] = common::create_sub_index_map(
1276 *x_index_map, sub_x_dofs, common::IndexMapOrder::any, true);
1277 sub_x_dof_index_map = std::make_shared<common::IndexMap>(std::move(map));
1278 subx_to_x_dofmap = std::move(new_to_old);
1279 }
1280
1281 // Create sub-geometry coordinates
1282 std::span<const T> x = geometry.x();
1283 std::int32_t sub_num_x_dofs = subx_to_x_dofmap.size();
1284 std::vector<T> sub_x(3 * sub_num_x_dofs);
1285 for (std::int32_t i = 0; i < sub_num_x_dofs; ++i)
1286 {
1287 std::copy_n(std::next(x.begin(), 3 * subx_to_x_dofmap[i]), 3,
1288 std::next(sub_x.begin(), 3 * i));
1289 }
1290
1291 // Create geometry to sub-geometry map
1292 std::vector<std::int32_t> x_to_subx_dof_map(
1293 x_index_map->size_local() + x_index_map->num_ghosts(), -1);
1294 for (std::size_t i = 0; i < subx_to_x_dofmap.size(); ++i)
1295 x_to_subx_dof_map[subx_to_x_dofmap[i]] = i;
1296
1297 // Create sub-geometry dofmap
1298 std::vector<std::int32_t> sub_x_dofmap;
1299 sub_x_dofmap.reserve(x_indices.size());
1300 std::ranges::transform(x_indices, std::back_inserter(sub_x_dofmap),
1301 [&x_to_subx_dof_map](auto x_dof)
1302 {
1303 assert(x_to_subx_dof_map[x_dof] != -1);
1304 return x_to_subx_dof_map[x_dof];
1305 });
1306
1307 // Sub-geometry coordinate element
1308 CellType sub_xcell = cell_entity_type(geometry.cmap().cell_shape(), dim, 0);
1309
1310 // Special handling of point meshes, as they only support constant
1311 // basis functions
1312 int degree = (sub_xcell == CellType::point) ? 0 : geometry.cmap().degree();
1313 fem::CoordinateElement<T> sub_cmap(sub_xcell, degree,
1314 geometry.cmap().variant());
1315
1316 // Sub-geometry input_global_indices
1317 const std::vector<std::int64_t>& igi = geometry.input_global_indices();
1318 std::vector<std::int64_t> sub_igi;
1319 sub_igi.reserve(subx_to_x_dofmap.size());
1320 std::ranges::transform(subx_to_x_dofmap, std::back_inserter(sub_igi),
1321 [&igi](auto sub_x_dof) { return igi[sub_x_dof]; });
1322
1323 // Create geometry
1324 return {Geometry(
1325 sub_x_dof_index_map,
1326 std::vector<std::vector<std::int32_t>>{std::move(sub_x_dofmap)},
1327 {sub_cmap}, std::move(sub_x), geometry.dim(), std::move(sub_igi)),
1328 std::move(subx_to_x_dofmap)};
1329}
1330
1340template <std::floating_point T>
1341std::tuple<Mesh<T>, std::vector<std::int32_t>, std::vector<std::int32_t>,
1342 std::vector<std::int32_t>>
1344 std::span<const std::int32_t> entities)
1345{
1346 // Create sub-topology
1347 mesh.topology_mutable()->create_connectivity(dim, 0);
1348 auto [topology, subentity_to_entity, subvertex_to_vertex]
1349 = mesh::create_subtopology(*mesh.topology(), dim, entities);
1350
1351 // Create sub-geometry
1352 const int tdim = mesh.topology()->dim();
1353 mesh.topology_mutable()->create_entities(dim);
1354 mesh.topology_mutable()->create_connectivity(dim, tdim);
1355 mesh.topology_mutable()->create_connectivity(tdim, dim);
1356 mesh.topology_mutable()->create_entity_permutations();
1357 auto [geometry, subx_to_x_dofmap]
1358 = mesh::create_subgeometry(mesh, dim, subentity_to_entity);
1359
1360 return {Mesh(mesh.comm(), std::make_shared<Topology>(std::move(topology)),
1361 std::move(geometry)),
1362 std::move(subentity_to_entity), std::move(subvertex_to_vertex),
1363 std::move(subx_to_x_dofmap)};
1364}
1365
1366} // namespace dolfinx::mesh
Definition CoordinateElement.h:38
ElementDofLayout create_dof_layout() const
Compute and return the dof layout.
Definition CoordinateElement.cpp:75
void permute_subentity_closure(std::span< std::int32_t > d, std::uint32_t cell_info, mesh::CellType entity_type, int entity_index) const
Given the closure DOFs of a cell sub-entity in reference ordering, this function computes the permut...
Definition CoordinateElement.cpp:64
Definition ElementDofLayout.h:30
const std::vector< std::vector< std::vector< int > > > & entity_closure_dofs_all() const
Definition ElementDofLayout.cpp:92
int num_entity_closure_dofs(int dim) const
Definition ElementDofLayout.cpp:68
Definition AdjacencyList.h:27
const std::vector< T > & array() const
Return contiguous array of links for all nodes (const version)
Definition AdjacencyList.h:128
const std::vector< std::int32_t > & offsets() const
Offset for each node in array() (const version)
Definition AdjacencyList.h:134
Geometry stores the geometry imposed on a mesh.
Definition Geometry.h:34
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:46
Requirements on function for geometry marking.
Definition utils.h:492
void reorder_list(std::span< T > list, std::span< const std::int32_t > nodemap)
Re-order the nodes of a fixed-degree adjacency list.
Definition utils.h:51
std::tuple< std::vector< std::int32_t >, std::vector< T >, std::vector< std::int32_t > > compute_vertex_coords_boundary(const mesh::Mesh< T > &mesh, int dim, std::span< const std::int32_t > facets)
Compute the coordinates of 'vertices' for entities of a given dimension that are attached to specifie...
Definition utils.h:82
std::pair< std::vector< T >, std::array< std::size_t, 2 > > compute_vertex_coords(const mesh::Mesh< T > &mesh)
The coordinates for all 'vertices' in the mesh.
Definition utils.h:446
std::vector< typename std::remove_reference_t< typename U::value_type > > distribute_data(MPI_Comm comm0, std::span< const std::int64_t > indices, MPI_Comm comm1, const U &x, int shape1)
Distribute rows of a rectangular data array to ranks where they are required (scalable version).
Definition MPI.h:680
int size(MPI_Comm comm)
Definition MPI.cpp:72
std::pair< IndexMap, std::vector< std::int32_t > > create_sub_index_map(const IndexMap &imap, std::span< const std::int32_t > indices, IndexMapOrder order=IndexMapOrder::any, bool allow_owner_change=false)
Create a new index map from a subset of indices in an existing index map.
Definition IndexMap.cpp:815
@ any
Allow arbitrary ordering of ghost indices in sub-maps.
Definition IndexMap.h:27
Finite element method functionality.
Definition assemble_expression_impl.h:23
Geometry data structures and algorithms.
Definition BoundingBoxTree.h:22
std::tuple< graph::AdjacencyList< std::int64_t >, std::vector< int >, std::vector< std::int64_t >, std::vector< int > > distribute(MPI_Comm comm, const graph::AdjacencyList< std::int64_t > &list, const graph::AdjacencyList< std::int32_t > &destinations)
Distribute adjacency list nodes to destination ranks.
Definition partition.cpp:38
Graph data structures and algorithms.
Definition dofmapbuilder.h:26
std::vector< std::int32_t > reorder_gps(const graph::AdjacencyList< std::int32_t > &graph)
Re-order a graph using the Gibbs-Poole-Stockmeyer algorithm.
Definition ordering.cpp:360
std::function< graph::AdjacencyList< std::int32_t >( MPI_Comm, int, const AdjacencyList< std::int64_t > &, bool)> partition_fn
Signature of functions for computing the parallel partitioning of a distributed graph.
Definition partition.h:31
AdjacencyList< std::int32_t > partition_graph(MPI_Comm comm, int nparts, const AdjacencyList< std::int64_t > &local_graph, bool ghosting)
Partition graph across processes using the default graph partitioner.
Definition partition.cpp:21
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
GhostMode
Enum for different partitioning ghost modes.
Definition utils.h:37
Topology create_topology(MPI_Comm comm, const std::vector< CellType > &cell_types, std::vector< std::span< const std::int64_t > > cells, std::vector< std::span< const std::int64_t > > original_cell_index, std::vector< std::span< const int > > ghost_owners, std::span< const std::int64_t > boundary_vertices)
Create a mesh topology.
Definition Topology.cpp:1017
Geometry< typename std::remove_reference_t< typename U::value_type > > create_geometry(const Topology &topology, const std::vector< fem::CoordinateElement< std::remove_reference_t< typename U::value_type > > > &elements, std::span< const std::int64_t > nodes, std::span< const std::int64_t > xdofs, const U &x, int dim, std::function< std::vector< int >(const graph::AdjacencyList< std::int32_t > &)> reorder_fn=nullptr)
Build Geometry from input data.
Definition Geometry.h:236
std::vector< T > cell_normals(const Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities)
Compute normal to given cell (viewed as embedded in 3D).
Definition utils.h:291
CellType cell_entity_type(CellType type, int d, int index)
Return type of cell for entity of dimension d at given entity index.
Definition cell_types.cpp:64
std::function< std::vector< std::int32_t >( const graph::AdjacencyList< std::int32_t > &)> CellReorderFunction
Function that reorders (locally) cells that are owned by this process. It takes the local mesh dual g...
Definition utils.h:214
std::tuple< Topology, std::vector< int32_t >, std::vector< int32_t > > create_subtopology(const Topology &topology, int dim, std::span< const std::int32_t > entities)
Create a topology for a subset of entities of a given topological dimension.
Definition Topology.cpp:1303
std::vector< std::int32_t > exterior_facet_indices(const Topology &topology, int facet_type_idx)
Compute the indices of all exterior facets that are owned by the caller.
Definition utils.cpp:58
std::vector< std::int32_t > locate_entities_boundary(const Mesh< T > &mesh, int dim, U marker)
Compute indices of all mesh entities that are attached to an owned boundary facet and evaluate to tru...
Definition utils.h:611
int num_cell_vertices(CellType type)
Definition cell_types.cpp:147
std::vector< T > h(const Mesh< T > &mesh, std::span< const std::int32_t > entities, int dim)
Compute greatest distance between any two vertices of the mesh entities (h).
Definition utils.h:239
std::pair< Geometry< T >, std::vector< int32_t > > create_subgeometry(const Mesh< T > &mesh, int dim, std::span< const std::int32_t > subentity_to_entity)
Create a sub-geometry from a mesh and a subset of mesh entities to be included.
Definition utils.h:1251
std::tuple< graph::AdjacencyList< std::int32_t >, std::vector< std::int64_t >, std::size_t, std::vector< std::int32_t > > build_local_dual_graph(std::span< const CellType > celltypes, const std::vector< std::span< const std::int64_t > > &cells)
Compute the local part of the dual graph (cell-cell connections via facets) and facets with only one ...
Definition graphbuild.cpp:354
std::pair< std::vector< std::int32_t >, std::array< std::size_t, 2 > > entities_to_geometry(const Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities, bool permute=false)
Compute the geometry degrees of freedom associated with the closure of a given set of cell entities.
Definition utils.h:687
std::vector< std::int32_t > compute_incident_entities(const Topology &topology, std::span< const std::int32_t > entities, int d0, int d1)
Compute incident entities.
Definition utils.cpp:123
std::vector< std::int64_t > extract_topology(CellType cell_type, const fem::ElementDofLayout &layout, std::span< const std::int64_t > cells)
Extract topology from cell data, i.e. extract cell vertices.
Definition utils.cpp:29
CellType
Cell type identifier.
Definition cell_types.h:22
std::vector< std::int32_t > locate_entities(const Mesh< T > &mesh, int dim, U marker, int entity_type_idx)
Compute indices of all mesh entities that evaluate to true for the provided geometric marking functio...
Definition utils.h:513
std::function< graph::AdjacencyList< std::int32_t >( MPI_Comm comm, int nparts, const std::vector< CellType > &cell_types, const std::vector< std::span< const std::int64_t > > &cells)> CellPartitionFunction
Signature for the cell partitioning function. Function that implement this interface compute the dest...
Definition utils.h:206
std::vector< T > compute_midpoints(const Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities)
Compute the midpoints for mesh entities of a given dimension.
Definition utils.h:408
Mesh< typename std::remove_reference_t< typename U::value_type > > create_mesh(MPI_Comm comm, MPI_Comm commt, std::vector< std::span< const std::int64_t > > cells, const std::vector< fem::CoordinateElement< typename std::remove_reference_t< typename U::value_type > > > &elements, MPI_Comm commg, const U &x, std::array< std::size_t, 2 > xshape, const CellPartitionFunction &partitioner, const CellReorderFunction &reorder_fn=graph::reorder_gps)
Create a distributed mesh::Mesh from mesh data and using the provided graph partitioning function for...
Definition utils.h:848
std::tuple< Mesh< T >, std::vector< std::int32_t >, std::vector< std::int32_t >, std::vector< std::int32_t > > create_submesh(const Mesh< T > &mesh, int dim, std::span< const std::int32_t > entities)
Create a new mesh consisting of a subset of entities in a mesh.
Definition utils.h:1343
CellPartitionFunction create_cell_partitioner(mesh::GhostMode ghost_mode=mesh::GhostMode::none, const graph::partition_fn &partfn=&graph::partition_graph)
Create a function that computes destination rank for mesh cells on this rank by applying the default ...
Definition utils.cpp:100
std::vector< std::int32_t > sort_by_perm(std::span< const T > x, std::size_t shape1)
Compute the permutation array that sorts a 2D array by row.
Definition sort.h:135
constexpr __radix_sort radix_sort
Radix sort.
Definition sort.h:122