Loading...
Searching...
No Matches
interval.h
// Copyright (C) 2024 Paul T. Kühner
//
// This file is part of DOLFINX (https://www.fenicsproject.org)
//
// SPDX-License-Identifier: LGPL-3.0-or-later
#pragma once
#include "dolfinx/mesh/Mesh.h"
#include "dolfinx/mesh/cell_types.h"
#include "dolfinx/mesh/utils.h"
#include "dolfinx/refinement/plaza.h"
#include <algorithm>
#include <concepts>
#include <cstddef>
#include <cstdint>
#include <mpi.h>
#include <optional>
#include <stdexcept>
#include <vector>
#include "dolfinx/refinement/option.h"
#include "dolfinx/refinement/utils.h"
namespace dolfinx::refinement::interval
{
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>> cells,
Option option)
{
bool compute_parent_facet = option_parent_facet(option);
bool compute_parent_cell = option_parent_cell(option);
if (compute_parent_facet)
throw std::runtime_error("Parent facet computation not yet supported!");
auto topology = mesh.topology();
assert(topology);
assert(topology->dim() == 1);
auto map_c = topology->index_map(1);
assert(map_c);
// TODO: creation of sharing ranks in external function? Also same
// code in use for plaza
// Get sharing ranks for each cell
graph::AdjacencyList<int> cell_ranks = map_c->index_to_dest_ranks();
// Create unique list of ranks that share cells (owners of ghosts plus
// ranks that ghost owned indices)
std::vector<int> ranks = cell_ranks.array();
std::ranges::sort(ranks);
auto to_remove = std::ranges::unique(ranks);
ranks.erase(to_remove.begin(), to_remove.end());
// Convert cell_ranks from global rank to to neighbourhood ranks
std::ranges::transform(cell_ranks.array(), cell_ranks.array().begin(),
[&ranks](auto r)
{
auto it = std::lower_bound(ranks.begin(),
ranks.end(), r);
assert(it != ranks.end() and *it == r);
return std::distance(ranks.begin(), it);
});
// Create refinement flag for cells
std::vector<std::int8_t> refinement_marker(
map_c->size_local() + map_c->num_ghosts(), !cells.has_value());
// Mark cells for refinement
std::vector<std::vector<std::int32_t>> marked_for_update(ranks.size());
if (cells)
{
std::ranges::for_each(
cells.value(),
[&refinement_marker, &cell_ranks, &marked_for_update](auto cell)
{
if (!refinement_marker[cell])
{
refinement_marker[cell] = true;
for (int rank : cell_ranks.links(cell))
marked_for_update[rank].push_back(cell);
}
});
}
// Create neighborhood communicator for vertex creation
MPI_Comm neighbor_comm;
MPI_Dist_graph_create_adjacent(
mesh.comm(), ranks.size(), ranks.data(), MPI_UNWEIGHTED, ranks.size(),
ranks.data(), MPI_UNWEIGHTED, MPI_INFO_NULL, false, &neighbor_comm);
// Communicate ghost cells that might have been marked. This is not
// necessary for a uniform refinement.
if (cells)
{
update_logical_edgefunction(neighbor_comm, marked_for_update,
refinement_marker, *map_c);
}
// Construct the new vertices
const auto [new_vertex_map, new_vertex_coords, xshape]
= create_new_vertices(neighbor_comm, cell_ranks, mesh, refinement_marker);
MPI_Comm_free(&neighbor_comm);
auto c_to_v = mesh.topology()->connectivity(1, 0);
assert(c_to_v);
// Get the count of cells to refine, note: we only consider non-ghost
// cells
const std::int32_t number_of_refined_cells
= std::count(refinement_marker.begin(),
std::next(refinement_marker.begin(),
mesh.topology()->index_map(1)->size_local()),
true);
// Produce local global indices, by padding out the previous index map
std::vector<std::int64_t> global_indices
= adjust_indices(*mesh.topology()->index_map(0), number_of_refined_cells);
// Build the topology on the new vertices
std::size_t refined_cell_count
= mesh.topology()->index_map(1)->size_local() + number_of_refined_cells;
std::vector<std::int64_t> cell_topology;
cell_topology.reserve(refined_cell_count * 2);
std::optional<std::vector<std::int32_t>> parent_cell(std::nullopt);
if (compute_parent_cell)
{
parent_cell.emplace();
parent_cell->reserve(refined_cell_count);
}
for (std::int32_t cell = 0; cell < map_c->size_local(); ++cell)
{
std::span vertices = c_to_v->links(cell);
assert(vertices.size() == 2);
// We consider a cell (defined by global vertices)
// a ----------- b
const std::int64_t a = global_indices[vertices[0]];
const std::int64_t b = global_indices[vertices[1]];
if (refinement_marker[cell])
{
// Find (global) index of new midpoint vertex:
// a --- c --- b
std::span nv = new_vertex_map.links(cell);
assert(nv.size() == 1);
const std::int64_t c = nv[0];
// Add new cells/edges to refined topology
cell_topology.insert(cell_topology.end(), {a, c, c, b});
if (compute_parent_cell)
parent_cell->insert(parent_cell->end(), {cell, cell});
}
else
{
// Copy the previous cell
cell_topology.insert(cell_topology.end(), {a, b});
if (compute_parent_cell)
parent_cell->push_back(cell);
}
}
assert(cell_topology.size() == 2 * refined_cell_count);
assert(!compute_parent_cell or parent_cell->size() == refined_cell_count);
std::vector<std::int32_t> offsets(refined_cell_count + 1);
std::ranges::generate(offsets, [i = 0]() mutable { return 2 * i++; });
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::nullopt};
}
} // namespace dolfinx::refinement::interval
utils.h
Functions supporting mesh operations.
dolfinx::fem::sparsitybuild::cells
void cells(la::SparsityPattern &pattern, std::array< std::span< const std::int32_t >, 2 > cells, std::array< std::reference_wrapper< const DofMap >, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition sparsitybuild.cpp:16
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
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