Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Remove inconsistent behaviour when building sparsity pattern from forms #2615

Merged
merged 3 commits into from
Apr 8, 2023
Merged

Remove inconsistent behaviour when building sparsity pattern from forms #2615

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
05d6284
Fix sparsity pattern for exterior facets
garth-wells Apr 7, 2023
81f9219
Simplify sparsity computation
garth-wells Apr 8, 2023
576c69d
Simplification
garth-wells Apr 8, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 19 additions & 74 deletions cpp/dolfinx/fem/petsc.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,91 +70,29 @@ Mat create_matrix_block(
bs_dofs[i].push_back(_V->dofmap()->bs());
}

std::shared_ptr mesh = V[0][0]->mesh();
assert(mesh);
auto topology = mesh->topology();
assert(topology);
const int tdim = topology->dim();

// Build sparsity pattern for each block
std::shared_ptr<const mesh::Mesh<T>> mesh;
std::vector<std::vector<std::unique_ptr<la::SparsityPattern>>> patterns(
V[0].size());
for (std::size_t row = 0; row < V[0].size(); ++row)
{
for (std::size_t col = 0; col < V[1].size(); ++col)
{
const std::array<std::shared_ptr<const common::IndexMap>, 2> index_maps
= {{V[0][row]->dofmap()->index_map, V[1][col]->dofmap()->index_map}};
const std::array bs = {V[0][row]->dofmap()->index_map_bs(),
V[1][col]->dofmap()->index_map_bs()};
if (const Form<PetscScalar, T>* form = a[row][col]; form)
{
// Create sparsity pattern for block
patterns[row].push_back(std::make_unique<la::SparsityPattern>(
mesh->comm(), index_maps, bs));

// Build sparsity pattern for block
assert(V[0][row]->dofmap());
assert(V[1][col]->dofmap());
std::array<std::reference_wrapper<const DofMap>, 2> dofmaps{
*V[0][row]->dofmap(), *V[1][col]->dofmap()};
assert(patterns[row].back());
auto& sp = patterns[row].back();
assert(sp);

if (form->num_integrals(IntegralType::cell) > 0)
{
auto map = topology->index_map(tdim);
assert(map);
std::vector<std::int32_t> c(map->size_local(), 0);
std::iota(c.begin(), c.end(), 0);
sparsitybuild::cells(*sp, c, dofmaps);
}

if (form->num_integrals(IntegralType::interior_facet) > 0)
{
// Loop over owned facets
mesh->topology_mutable()->create_entities(tdim - 1);
auto f_to_c = topology->connectivity(tdim - 1, tdim);
if (!f_to_c)
{
throw std::runtime_error(
"Facet-cell connectivity has not been computed.");
}
auto map = topology->index_map(tdim - 1);
assert(map);
std::vector<std::int32_t> facets;
facets.reserve(2 * map->size_local());
for (int f = 0; f < map->size_local(); ++f)
if (auto cells = f_to_c->links(f); cells.size() == 2)
facets.insert(facets.end(), {cells[0], cells[1]});
sparsitybuild::interior_facets(*sp, facets, dofmaps);
}

if (form->num_integrals(IntegralType::exterior_facet) > 0)
{
// Loop over owned facets
mesh->topology_mutable()->create_entities(tdim - 1);
auto connectivity = topology->connectivity(tdim - 1, tdim);
if (!connectivity)
{
throw std::runtime_error(
"Facet-cell connectivity has not been computed.");
}
auto map = topology->index_map(tdim - 1);
assert(map);
std::vector<std::int32_t> cells;
for (int f = 0; f < map->size_local(); ++f)
if (auto c = connectivity->links(f); c.size() == 1)
cells.push_back(c[0]);
sparsitybuild::cells(*sp, cells, dofmaps);
}
create_sparsity_pattern(*form)));
if (!mesh)
mesh = form->mesh();
}
else
patterns[row].push_back(nullptr);
}
}

if (!mesh)
throw std::runtime_error("Could not find a Mesh.");

// Compute offsets for the fields
std::array<std::vector<std::pair<
std::reference_wrapper<const common::IndexMap>, int>>,
Expand Down Expand Up @@ -240,9 +178,9 @@ Mat create_matrix_block(
return A;
}

/// Create nested (MatNest) matrix
/// @brief Create nested (MatNest) matrix.
///
/// The caller is responsible for destroying the Mat object
/// @note The caller is responsible for destroying the Mat object.
template <std::floating_point T>
Mat create_matrix_nest(
const std::vector<std::vector<const Form<PetscScalar, T>*>>& a,
Expand All @@ -257,21 +195,28 @@ Mat create_matrix_nest(
_types = types;

// Loop over each form and create matrix
const int rows = a.size();
const int cols = a.front().size();
int rows = a.size();
int cols = a.front().size();
std::vector<Mat> mats(rows * cols, nullptr);
std::shared_ptr<const mesh::Mesh<T>> mesh;
for (int i = 0; i < rows; ++i)
{
for (int j = 0; j < cols; ++j)
{
if (const Form<PetscScalar, T>* form = a[i][j]; form)
{
mats[i * cols + j] = create_matrix(*form, _types[i][j]);
mesh = form->mesh();
}
}
}

if (!mesh)
throw std::runtime_error("Could not find a Mesh.");

// Initialise block (MatNest) matrix
Mat A;
MatCreate(V[0][0]->mesh()->comm(), &A);
MatCreate(mesh->comm(), &A);
MatSetType(A, MATNEST);
MatNestSetSubMats(A, rows, nullptr, cols, nullptr, mats.data());
MatSetUp(A);
Expand Down