83 std::span<const std::int32_t> facets)
85 auto topology =
mesh.topology();
87 const int tdim = topology->dim();
90 throw std::runtime_error(
91 "Cannot use mesh::locate_entities_boundary (boundary) for cells.");
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;
99 auto f_to_v = topology->connectivity(tdim - 1, 0);
101 auto f_to_e = topology->connectivity(tdim - 1, dim);
103 for (
auto f : facets)
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());
113 std::ranges::sort(vertices);
114 auto [unique_end, range_end] = std::ranges::unique(vertices);
115 vertices.erase(unique_end, range_end);
119 std::ranges::sort(entities);
120 auto [unique_end, range_end] = std::ranges::unique(entities);
121 entities.erase(unique_end, range_end);
126 auto x_dofmap =
mesh.geometry().dofmap();
127 std::span<const T> x_nodes =
mesh.geometry().x();
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);
134 auto c_to_v = topology->connectivity(tdim, 0);
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)
140 const std::int32_t v = vertices[i];
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);
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;
155 return {std::move(entities), std::move(x_vertices), std::move(vertex_to_pos)};
292 std::span<const std::int32_t> entities)
294 if (entities.empty())
295 return std::vector<T>();
297 auto topology =
mesh.topology();
299 if (topology->cell_type() == CellType::prism and dim == 2)
301 throw std::runtime_error(
302 "Cell normal computation for prism cells not yet supported.");
305 const int gdim =
mesh.geometry().dim();
309 std::span<const T> x =
mesh.geometry().x();
310 const auto [geometry_entities, eshape]
313 std::vector<T> n(entities.size() * 3);
316 case CellType::interval:
319 throw std::invalid_argument(
"Interval cell normal undefined in 3D.");
320 for (std::size_t i = 0; i < entities.size(); ++i)
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)};
330 std::ranges::transform(p[1], p[0], t.begin(),
331 [](
auto x,
auto y) { return x - y; });
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;
341 case CellType::triangle:
343 for (std::size_t i = 0; i < entities.size(); ++i)
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)};
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; });
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; });
369 case CellType::quadrilateral:
372 for (std::size_t i = 0; i < entities.size(); ++i)
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)};
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; });
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; });
399 throw std::invalid_argument(
400 "cell_normal not supported for this cell type.");
688 std::span<const std::int32_t> entities,
689 bool permute =
false)
691 auto topology =
mesh.topology();
693 CellType cell_type = topology->cell_type();
694 if (cell_type == CellType::prism and dim == 2)
696 throw std::runtime_error(
697 "mesh::entities_to_geometry for prism cells not yet supported.");
700 const int tdim = topology->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};
713 const std::vector<std::vector<std::vector<int>>>& closure_dofs_all
719 for (std::int32_t c : entities)
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]);
727 return {std::move(entity_xdofs), eshape};
732 auto e_to_c = topology->connectivity(dim, tdim);
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));
740 auto c_to_e = topology->connectivity(tdim, dim);
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));
749 std::span<const std::uint32_t> cell_info;
751 cell_info = std::span(
mesh.topology()->get_cell_permutation_info());
753 for (std::int32_t e : entities)
756 assert(!e_to_c->links(e).empty());
757 std::int32_t c = e_to_c->links(e).front();
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);
767 std::vector<std::int32_t> closure_dofs(closure_dofs_all[dim][local_entity]);
773 entity_type, local_entity);
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]);
782 return {std::move(entity_xdofs), eshape};
849 MPI_Comm comm, MPI_Comm commt,
850 std::vector<std::span<const std::int64_t>> cells,
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,
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(); });
872 std::int32_t num_cell_types = cells.size();
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);
880 spdlog::info(
"Using partitioner with cell data ({} cell types)",
883 if (commt != MPI_COMM_NULL)
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)
891 tspan[i] = std::span(t[i]);
893 dest = partitioner(commt, size, celltypes, tspan);
896 std::int32_t cell_offset = 0;
897 for (std::int32_t i = 0; i < num_cell_types; ++i)
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;
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; });
914 cell_offset += num_cells;
918 std::vector<int> src_ranks;
919 std::tie(cells1[i], src_ranks, original_idx1[i], ghost_owners[i])
921 {num_cells, num_cell_nodes}, dest_i);
922 spdlog::debug(
"Got {} cells from distribution", cells1[i].size());
928 std::int64_t num_owned = 0;
929 for (std::int32_t i = 0; i < num_cell_types; ++i)
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();
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)
943 std::iota(original_idx1[i].begin(), original_idx1[i].end(),
945 global_offset += original_idx1[i].size();
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)
955 spdlog::info(
"Extract basic topology: {}->{}", cells1[i].size(),
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,
970 spdlog::info(
"Build local dual graphs, re-order cells, and compute process "
971 "boundary vertices.");
973 std::vector<std::pair<std::vector<std::int64_t>,
int>> facets;
976 std::vector<std::span<const std::int64_t>> cells1_v_local;
977 for (std::size_t i = 0; i < celltypes.size(); ++i)
980 std::size_t num_owned_cells
982 cells1_v_local.emplace_back(cells1_v[i].data(),
986 auto [
graph, unmatched_facets, max_v, _facet_attached_cells]
988 std::vector{cells1_v_local.back()});
991 facets.emplace_back(std::move(unmatched_facets), max_v);
994 const std::vector<std::int32_t> remap = reorder_fn(
graph);
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());
1002 for (std::size_t j = 0; j < remap.size(); ++j)
1003 _original_idx[remap[j]] = orig_idx[j];
1005 original_idx1[i] = _original_idx;
1012 std::span(cells1[i].data(), remap.size() * doflayouts[i].num_dofs()),
1016 if (facets.size() == 1)
1018 std::vector<std::int64_t>& vertices = facets.front().first;
1021 std::ranges::sort(vertices);
1022 auto [unique_end, range_end] = std::ranges::unique(vertices);
1023 vertices.erase(unique_end, range_end);
1028 if (!vertices.empty() and vertices.front() == -1)
1029 vertices.erase(vertices.begin());
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; })
1044 for (
const auto& [v_data, num_v] : facets)
1046 for (
auto it = v_data.begin(); it != v_data.end(); it += num_v)
1048 facets0.insert(facets0.end(), it, std::next(it, num_v));
1049 facets0.insert(facets0.end(), max_v - num_v, -1);
1055 std::span<const std::int64_t>(facets0), max_v);
1059 std::vector<std::int64_t> vertices;
1060 auto it = perm.begin();
1061 while (it != perm.end())
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()));
1070 auto it1 = std::find_if_not(
1072 [f, max_v, it0 = facets0.begin()](
auto p) ->
bool
1074 return std::equal(f.begin(), f.end(), std::next(it0, p * max_v));
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.");
1088 std::ranges::sort(vertices);
1089 auto [unique_end, range_end] = std::ranges::unique(vertices);
1090 vertices.erase(unique_end, range_end);
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);
1100 spdlog::debug(
"Got {} boundary vertices", boundary_v.size());
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); });
1114 ghost_owners_span, boundary_v);
1118 for (
int i = 0; i < num_cell_types; ++i)
1120 for (
int e = 1; e < topology.dim(); ++e)
1122 if (doflayouts[i].num_entity_dofs(e) > 0)
1123 topology.create_entities(e);
1126 if (elements[i].needs_dof_permutations())
1127 topology.create_entity_permutations();
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());
1139 auto [unique_end, range_end] = std::ranges::unique(nodes1);
1140 nodes1.erase(unique_end, range_end);
1147 =
create_geometry(topology, elements, nodes1, nodes2, coords, xshape[1]);
1149 return Mesh(comm, std::make_shared<Topology>(std::move(topology)),
1252 std::span<const std::int32_t> subentity_to_entity)
1260 const std::vector<std::int32_t> x_indices
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);
1269 auto x_index_map =
geometry.index_map();
1270 assert(x_index_map);
1272 std::shared_ptr<common::IndexMap> sub_x_dof_index_map;
1273 std::vector<std::int32_t> subx_to_x_dofmap;
1277 sub_x_dof_index_map = std::make_shared<common::IndexMap>(std::move(map));
1278 subx_to_x_dofmap = std::move(new_to_old);
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)
1287 std::copy_n(std::next(x.begin(), 3 * subx_to_x_dofmap[i]), 3,
1288 std::next(sub_x.begin(), 3 * i));
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;
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)
1303 assert(x_to_subx_dof_map[x_dof] != -1);
1304 return x_to_subx_dof_map[x_dof];
1312 int degree = (sub_xcell == CellType::point) ? 0 :
geometry.cmap().degree();
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]; });
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)};