Loading...
Searching...
No Matches
plaza.h
// Copyright (C) 2014-2018 Chris Richardson
//
// This file is part of DOLFINx (https://www.fenicsproject.org)
//
// SPDX-License-Identifier: LGPL-3.0-or-later
#include "dolfinx/graph/AdjacencyList.h"
#include "dolfinx/mesh/Mesh.h"
#include "dolfinx/mesh/Topology.h"
#include "dolfinx/mesh/utils.h"
#include "option.h"
#include "utils.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <optional>
#include <span>
#include <tuple>
#include <utility>
#include <vector>
#pragma once
namespace dolfinx::refinement::plaza
{
namespace impl
{
template <int tdim>
auto compute_parent_facets(std::span<const std::int32_t> simplex_set)
{
static_assert(tdim == 2 or tdim == 3);
assert(simplex_set.size() % (tdim + 1) == 0);
using parent_facet_t
= std::conditional_t<tdim == 2, std::array<std::int8_t, 12>,
std::array<std::int8_t, 32>>;
parent_facet_t parent_facet;
parent_facet.fill(-1);
assert(simplex_set.size() <= parent_facet.size());
// Index lookups in 'indices' for the child vertices that occur on
// each parent facet in 2D and 3D. In 2D each edge has 3 child
// vertices, and in 3D each triangular facet has six child vertices.
constexpr std::array<std::array<int, 3>, 3> facet_table_2d{
{{1, 2, 3}, {0, 2, 4}, {0, 1, 5}}};
constexpr std::array<std::array<int, 6>, 4> facet_table_3d{
{{1, 2, 3, 4, 5, 6},
{0, 2, 3, 4, 7, 8},
{0, 1, 3, 5, 7, 9},
{0, 1, 2, 6, 8, 9}}};
const int ncells = simplex_set.size() / (tdim + 1);
for (int fpi = 0; fpi < (tdim + 1); ++fpi)
{
// For each child cell, consider all facets
for (int cc = 0; cc < ncells; ++cc)
{
for (int fci = 0; fci < (tdim + 1); ++fci)
{
// Indices of all vertices on child facet, sorted
std::array<int, tdim> cf, set_output;
int num_common_vertices;
if constexpr (tdim == 2)
{
for (int j = 0; j < tdim; ++j)
cf[j] = simplex_set[cc * 3 + facet_table_2d[fci][j]];
std::ranges::sort(cf);
auto [last1, last2, it_last] = std::ranges::set_intersection(
facet_table_2d[fpi], cf, set_output.begin());
num_common_vertices = std::distance(set_output.begin(), it_last);
}
else
{
for (int j = 0; j < tdim; ++j)
cf[j] = simplex_set[cc * 4 + facet_table_3d[fci][j]];
std::ranges::sort(cf);
auto [last1, last2, it_last] = std::ranges::set_intersection(
facet_table_3d[fpi], cf, set_output.begin());
num_common_vertices = std::distance(set_output.begin(), it_last);
}
if (num_common_vertices == tdim)
{
assert(parent_facet[cc * (tdim + 1) + fci] == -1);
// Child facet "fci" of cell cc, lies on parent facet "fpi"
parent_facet[cc * (tdim + 1) + fci] = fpi;
}
}
}
}
return parent_facet;
}
std::pair<std::array<std::int32_t, 32>, std::size_t>
get_simplices(std::span<const std::int64_t> indices,
std::span<const std::int32_t> longest_edge, int tdim,
bool uniform);
void enforce_rules(MPI_Comm comm, const graph::AdjacencyList<int>& shared_edges,
std::span<std::int8_t> marked_edges,
const mesh::Topology& topology,
std::span<const std::int32_t> long_edge);
template <std::floating_point T>
std::pair<std::vector<std::int32_t>, std::vector<std::int8_t>>
face_long_edge(const mesh::Mesh<T>& mesh)
{
const int tdim = mesh.topology()->dim();
// FIXME: cleanup these calls? Some of the happen internally again.
mesh.topology_mutable()->create_entities(1);
mesh.topology_mutable()->create_entities(2);
mesh.topology_mutable()->create_connectivity(2, 1);
mesh.topology_mutable()->create_connectivity(1, tdim);
mesh.topology_mutable()->create_connectivity(tdim, 2);
std::int64_t num_faces = mesh.topology()->index_map(2)->size_local()
+ mesh.topology()->index_map(2)->num_ghosts();
// Storage for face-local index of longest edge
std::vector<std::int32_t> long_edge(num_faces);
std::vector<std::int8_t> edge_ratio_ok;
// Check mesh face quality (may be used in 2D to switch to "uniform"
// refinement)
const T min_ratio = std::sqrt(2.0) / 2.0;
if (tdim == 2)
edge_ratio_ok.resize(num_faces);
auto x_dofmap = mesh.geometry().dofmap();
auto c_to_v = mesh.topology()->connectivity(tdim, 0);
assert(c_to_v);
auto e_to_c = mesh.topology()->connectivity(1, tdim);
assert(e_to_c);
auto e_to_v = mesh.topology()->connectivity(1, 0);
assert(e_to_v);
// Store all edge lengths in Mesh to save recalculating for each Face
auto map_e = mesh.topology()->index_map(1);
assert(map_e);
std::vector<T> edge_length(map_e->size_local() + map_e->num_ghosts());
for (std::size_t e = 0; e < edge_length.size(); ++e)
{
// Get first attached cell
auto cells = e_to_c->links(e);
assert(!cells.empty());
auto cell_vertices = c_to_v->links(cells.front());
auto edge_vertices = e_to_v->links(e);
// Find local index of edge vertices in the cell geometry map
auto it0 = std::find(cell_vertices.begin(), cell_vertices.end(),
edge_vertices[0]);
assert(it0 != cell_vertices.end());
const std::size_t local0 = std::distance(cell_vertices.begin(), it0);
auto it1 = std::find(cell_vertices.begin(), cell_vertices.end(),
edge_vertices[1]);
assert(it1 != cell_vertices.end());
const std::size_t local1 = std::distance(cell_vertices.begin(), it1);
auto x_dofs = md::submdspan(x_dofmap, cells.front(), md::full_extent);
std::span<const T, 3> x0(mesh.geometry().x().data() + 3 * x_dofs[local0],
3);
std::span<const T, 3> x1(mesh.geometry().x().data() + 3 * x_dofs[local1],
3);
// Compute length of edge between vertex x0 and x1
edge_length[e] = std::sqrt(std::transform_reduce(
x0.begin(), x0.end(), x1.begin(), 0.0, std::plus<>(),
[](auto x0, auto x1) { return (x0 - x1) * (x0 - x1); }));
}
// Get longest edge of each face
auto f_to_v = mesh.topology()->connectivity(2, 0);
assert(f_to_v);
auto f_to_e = mesh.topology()->connectivity(2, 1);
assert(f_to_e);
const std::vector global_indices
= mesh.topology()->index_map(0)->global_indices();
for (int f = 0; f < f_to_v->num_nodes(); ++f)
{
auto face_edges = f_to_e->links(f);
std::int32_t imax = 0;
T max_len = 0.0;
T min_len = std::numeric_limits<T>::max();
for (int i = 0; i < 3; ++i)
{
const T e_len = edge_length[face_edges[i]];
min_len = std::min(e_len, min_len);
if (e_len > max_len)
{
max_len = e_len;
imax = i;
}
else if (tdim == 3 and e_len == max_len)
{
// If edges are the same length, compare global index of
// opposite vertex. Only important so that tetrahedral faces
// have a matching refinement pattern across processes.
auto vertices = f_to_v->links(f);
const int vmax = vertices[imax];
const int vi = vertices[i];
if (global_indices[vi] > global_indices[vmax])
imax = i;
}
}
// Only save edge ratio in 2D
if (tdim == 2)
edge_ratio_ok[f] = (min_len / max_len >= min_ratio);
long_edge[f] = face_edges[imax];
}
return std::pair(std::move(long_edge), std::move(edge_ratio_ok));
}
template <std::floating_point T>
std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
std::optional<std::vector<std::int8_t>>>
compute_refinement(MPI_Comm neighbor_comm,
std::span<const std::int8_t> marked_edges,
const graph::AdjacencyList<int>& shared_edges,
const mesh::Mesh<T>& mesh,
std::span<const std::int32_t> long_edge,
std::span<const std::int8_t> edge_ratio_ok, Option option)
{
int tdim = mesh.topology()->dim();
int num_cell_edges = tdim * 3 - 3;
int num_cell_vertices = tdim + 1;
bool compute_facets = option_parent_facet(option);
bool compute_parent_cell = option_parent_cell(option);
// Make new vertices in parallel
const auto [new_vertex_map, new_vertex_coords, xshape]
= create_new_vertices(neighbor_comm, shared_edges, mesh, marked_edges);
std::optional<std::vector<std::int32_t>> parent_cell(std::nullopt);
if (compute_parent_cell)
parent_cell.emplace();
std::optional<std::vector<std::int8_t>> parent_facet(std::nullopt);
if (compute_facets)
parent_facet.emplace();
std::vector<std::int64_t> indices(num_cell_vertices + num_cell_edges);
std::vector<std::int32_t> simplex_set;
auto map_c = mesh.topology()->index_map(tdim);
assert(map_c);
auto c_to_v = mesh.topology()->connectivity(tdim, 0);
assert(c_to_v);
auto c_to_e = mesh.topology()->connectivity(tdim, 1);
assert(c_to_e);
auto c_to_f = mesh.topology()->connectivity(tdim, 2);
assert(c_to_f);
std::int32_t num_new_vertices_local = std::count(
marked_edges.begin(),
marked_edges.begin() + mesh.topology()->index_map(1)->size_local(), true);
std::vector<std::int64_t> global_indices
= adjust_indices(*mesh.topology()->index_map(0), num_new_vertices_local);
const std::int32_t num_cells = map_c->size_local();
// Iterate over all cells, and refine if cell has a marked edge
std::vector<std::int64_t> cell_topology;
for (int c = 0; c < num_cells; ++c)
{
// Create vector of indices in the order [vertices][edges], 3+3 in
// 2D, 4+6 in 3D
// Copy vertices
auto vertices = c_to_v->links(c);
for (std::size_t v = 0; v < vertices.size(); ++v)
indices[v] = global_indices[vertices[v]];
// Get cell-local indices of marked edges
auto edges = c_to_e->links(c);
bool no_edge_marked = true;
for (std::size_t ei = 0; ei < edges.size(); ++ei)
{
if (marked_edges[edges[ei]])
{
no_edge_marked = false;
auto nv = new_vertex_map.links(edges[ei]);
assert(nv.size() == 1);
indices[num_cell_vertices + ei] = nv[0];
}
else
indices[num_cell_vertices + ei] = -1;
}
if (no_edge_marked)
{
// Copy over existing cell to new topology
for (auto v : vertices)
cell_topology.push_back(global_indices[v]);
if (compute_parent_cell)
parent_cell->push_back(c);
if (compute_facets)
{
if (tdim == 3)
parent_facet->insert(parent_facet->end(), {0, 1, 2, 3});
else
parent_facet->insert(parent_facet->end(), {0, 1, 2});
}
}
else
{
// Need longest edges of each face in cell local indexing. NB in
// 2D the face is the cell itself, and there is just one entry.
std::vector<std::int32_t> longest_edge;
for (auto f : c_to_f->links(c))
longest_edge.push_back(long_edge[f]);
// Convert to cell local index
for (std::int32_t& p : longest_edge)
{
for (std::size_t ej = 0; ej < edges.size(); ++ej)
{
if (p == edges[ej])
{
p = ej;
break;
}
}
}
const bool uniform = (tdim == 2) ? edge_ratio_ok[c] : false;
const auto [simplex_set_b, simplex_set_size]
= get_simplices(indices, longest_edge, tdim, uniform);
std::span<const std::int32_t> simplex_set(simplex_set_b.data(),
simplex_set_size);
// Save parent index
const std::int32_t ncells = simplex_set.size() / num_cell_vertices;
if (compute_parent_cell)
{
for (std::int32_t i = 0; i < ncells; ++i)
parent_cell->push_back(c);
}
if (compute_facets)
{
if (tdim == 3)
{
auto npf = compute_parent_facets<3>(simplex_set);
parent_facet->insert(parent_facet->end(), npf.begin(),
std::next(npf.begin(), simplex_set.size()));
}
else
{
auto npf = compute_parent_facets<2>(simplex_set);
parent_facet->insert(parent_facet->end(), npf.begin(),
std::next(npf.begin(), simplex_set.size()));
}
}
// Convert from cell local index to mesh index and add to cells
for (std::int32_t v : simplex_set)
cell_topology.push_back(indices[v]);
}
}
assert(cell_topology.size() % num_cell_vertices == 0);
std::vector<std::int32_t> offsets(
cell_topology.size() / num_cell_vertices + 1, 0);
for (std::size_t i = 0; i < offsets.size() - 1; ++i)
offsets[i + 1] = offsets[i] + num_cell_vertices;
graph::AdjacencyList cell_adj(std::move(cell_topology), std::move(offsets));
return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
std::move(parent_cell), std::move(parent_facet)};
}
} // namespace impl
template <std::floating_point T>
std::tuple<graph::AdjacencyList<std::int64_t>, std::vector<T>,
std::array<std::size_t, 2>, std::optional<std::vector<std::int32_t>>,
std::optional<std::vector<std::int8_t>>>
compute_refinement_data(const mesh::Mesh<T>& mesh,
std::optional<std::span<const std::int32_t>> edges,
Option option)
{
common::Timer t0("PLAZA: refine");
auto topology = mesh.topology();
assert(topology);
if (topology->cell_type() != mesh::CellType::triangle
and topology->cell_type() != mesh::CellType::tetrahedron)
{
throw std::runtime_error("Cell type not supported");
}
auto map_e = topology->index_map(1);
if (!map_e)
throw std::runtime_error("Edges must be initialised");
// Get sharing ranks for each edge
graph::AdjacencyList<int> edge_ranks = map_e->index_to_dest_ranks();
// Create unique list of ranks that share edges (owners of ghosts plus
// ranks that ghost owned indices)
std::vector<int> ranks(edge_ranks.array().begin(), edge_ranks.array().end());
std::ranges::sort(ranks);
auto [unique_end, range_end] = std::ranges::unique(ranks);
ranks.erase(unique_end, range_end);
// Convert edge_ranks from global rank to to neighbourhood ranks
std::ranges::transform(edge_ranks.array(), edge_ranks.array().begin(),
[&ranks](auto r)
{
auto it = std::ranges::lower_bound(ranks, r);
assert(it != ranks.end() and *it == r);
return std::distance(ranks.begin(), it);
});
// Get number of neighbors
std::vector<std::int8_t> marked_edges(
map_e->size_local() + map_e->num_ghosts(), !edges.has_value());
std::vector<std::vector<std::int32_t>> marked_for_update(ranks.size());
if (edges)
{
for (auto edge : edges.value())
{
if (!marked_edges[edge])
{
marked_edges[edge] = true;
// If it is a shared edge, add all sharing neighbors to update set
for (int rank : edge_ranks.links(edge))
marked_for_update[rank].push_back(edge);
}
}
}
MPI_Comm comm;
MPI_Dist_graph_create_adjacent(mesh.comm(), ranks.size(), ranks.data(),
MPI_UNWEIGHTED, ranks.size(), ranks.data(),
MPI_UNWEIGHTED, MPI_INFO_NULL, false, &comm);
// Communicate any shared edges
update_logical_edgefunction(comm, marked_for_update, marked_edges, *map_e);
// Enforce rules about refinement (i.e. if any edge is marked in a
// triangle, then the longest edge must also be marked).
const auto [long_edge, edge_ratio_ok] = impl::face_long_edge(mesh);
impl::enforce_rules(comm, edge_ranks, marked_edges, *topology, long_edge);
auto [cell_adj, new_vertex_coords, xshape, parent_cell, parent_facet]
= impl::compute_refinement(comm, marked_edges, edge_ranks, mesh,
long_edge, edge_ratio_ok, option);
MPI_Comm_free(&comm);
return {std::move(cell_adj), std::move(new_vertex_coords), xshape,
std::move(parent_cell), std::move(parent_facet)};
}
} // namespace dolfinx::refinement::plaza
dolfinx::common::Timer
Timer for measuring and logging elapsed time durations.
Definition Timer.h:40
dolfinx::graph::AdjacencyList
This class provides a static adjacency list data structure.
Definition AdjacencyList.h:38
dolfinx::graph::AdjacencyList::array
const std::vector< LinkData > & array() const
Return contiguous array of links for all nodes (const version).
Definition AdjacencyList.h:173
dolfinx::graph::AdjacencyList::links
std::span< LinkData > links(std::size_t node)
Get the links (edges) for given node.
Definition AdjacencyList.h:156
dolfinx::mesh::Mesh
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition Mesh.h:23
dolfinx::mesh::Topology
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition Topology.h:41
utils.h
Functions supporting mesh operations.
dolfinx::mesh
Mesh data structures and algorithms on meshes.
Definition DofMap.h:32
dolfinx::refinement::plaza
Plaza mesh refinement.
Definition plaza.h:30
dolfinx::refinement::plaza::compute_refinement_data
std::tuple< graph::AdjacencyList< std::int64_t >, std::vector< T >, std::array< std::size_t, 2 >, std::optional< std::vector< std::int32_t > >, std::optional< std::vector< std::int8_t > > > compute_refinement_data(const mesh::Mesh< T > &mesh, std::optional< std::span< const std::int32_t > > edges, Option option)
Definition plaza.h:462
dolfinx::refinement::create_new_vertices
std::tuple< graph::AdjacencyList< std::int64_t >, std::vector< T >, std::array< std::size_t, 2 > > create_new_vertices(MPI_Comm comm, const graph::AdjacencyList< int > &shared_edges, const mesh::Mesh< T > &mesh, std::span< const std::int8_t > marked_edges)
Add new vertex for each marked edge, and create new_vertex_coordinates and global_edge->new_vertex ma...
Definition utils.h:148
dolfinx::refinement::option_parent_cell
constexpr bool option_parent_cell(Option a)
Check if parent_cell flag is set.
Definition option.h:45
dolfinx::refinement::Option
Option
Options for data to compute during mesh refinement.
Definition option.h:16
dolfinx::refinement::Option::parent_cell
@ parent_cell
Definition option.h:20
dolfinx::refinement::Option::parent_facet
@ parent_facet
Definition option.h:18
dolfinx::refinement::option_parent_facet
constexpr bool option_parent_facet(Option a)
Check if parent_facet flag is set.
Definition option.h:36
dolfinx::refinement::adjust_indices
std::vector< std::int64_t > adjust_indices(const common::IndexMap &map, std::int32_t n)
Given an index map, add "n" extra indices at the end of local range.
Definition utils.cpp:103
dolfinx::refinement::update_logical_edgefunction
void update_logical_edgefunction(MPI_Comm comm, const std::vector< std::vector< std::int32_t > > &marked_for_update, std::span< std::int8_t > marked_edges, const common::IndexMap &map)
Communicate edge markers between processes that share edges.
Definition utils.cpp:47