From 05d6284f4dda077d469b056fc90dfb685803423d Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Fri, 7 Apr 2023 10:50:50 +0200 Subject: [PATCH 1/3] Fix sparsity pattern for exterior facets --- cpp/dolfinx/fem/petsc.h | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h index 92b41e4b563..cefc409042a 100644 --- a/cpp/dolfinx/fem/petsc.h +++ b/cpp/dolfinx/fem/petsc.h @@ -133,20 +133,21 @@ Mat create_matrix_block( if (form->num_integrals(IntegralType::exterior_facet) > 0) { - // Loop over owned facets mesh->topology_mutable()->create_entities(tdim - 1); + const std::vector facets = mesh::exterior_facet_indices(*topology); 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 cells; - for (int f = 0; f < map->size_local(); ++f) - if (auto c = connectivity->links(f); c.size() == 1) - cells.push_back(c[0]); + for (auto f : facets) + { + auto c = connectivity->links(f); + assert(c.size() == 1); + cells.push_back(c[0]); + } sparsitybuild::cells(*sp, cells, dofmaps); } } From 81f9219817280dc38d58ea570c121ef3fdb79a9d Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 8 Apr 2023 10:11:29 +0200 Subject: [PATCH 2/3] Simplify sparsity computation --- cpp/dolfinx/fem/petsc.h | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h index cefc409042a..18cbd154b26 100644 --- a/cpp/dolfinx/fem/petsc.h +++ b/cpp/dolfinx/fem/petsc.h @@ -72,9 +72,9 @@ Mat create_matrix_block( std::shared_ptr mesh = V[0][0]->mesh(); assert(mesh); - auto topology = mesh->topology(); - assert(topology); - const int tdim = topology->dim(); + // auto topology = mesh->topology(); + // assert(topology); + // const int tdim = topology->dim(); // Build sparsity pattern for each block std::vector>> patterns( @@ -83,16 +83,22 @@ Mat create_matrix_block( { for (std::size_t col = 0; col < V[1].size(); ++col) { - const std::array, 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()}; + // const std::array, 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* form = a[row][col]; form) { + // Create sparsity pattern for block + patterns[row].push_back(std::make_unique( + create_sparsity_pattern(*form))); + + /* // Create sparsity pattern for block patterns[row].push_back(std::make_unique( mesh->comm(), index_maps, bs)); + // Build sparsity pattern for block assert(V[0][row]->dofmap()); assert(V[1][col]->dofmap()); @@ -150,6 +156,7 @@ Mat create_matrix_block( } sparsitybuild::cells(*sp, cells, dofmaps); } + */ } else patterns[row].push_back(nullptr); From 576c69d469c34c2d8bac0b3b38a07f735c644df0 Mon Sep 17 00:00:00 2001 From: "Garth N. Wells" Date: Sat, 8 Apr 2023 10:29:28 +0200 Subject: [PATCH 3/3] Simplification --- cpp/dolfinx/fem/petsc.h | 99 ++++++++--------------------------------- 1 file changed, 18 insertions(+), 81 deletions(-) diff --git a/cpp/dolfinx/fem/petsc.h b/cpp/dolfinx/fem/petsc.h index 18cbd154b26..2d5ad8319bd 100644 --- a/cpp/dolfinx/fem/petsc.h +++ b/cpp/dolfinx/fem/petsc.h @@ -70,99 +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> mesh; std::vector>> 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, 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* form = a[row][col]; form) { - // Create sparsity pattern for block patterns[row].push_back(std::make_unique( create_sparsity_pattern(*form))); - - /* - // Create sparsity pattern for block - patterns[row].push_back(std::make_unique( - mesh->comm(), index_maps, bs)); - - - // Build sparsity pattern for block - assert(V[0][row]->dofmap()); - assert(V[1][col]->dofmap()); - std::array, 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 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 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) - { - mesh->topology_mutable()->create_entities(tdim - 1); - const std::vector facets = mesh::exterior_facet_indices(*topology); - auto connectivity = topology->connectivity(tdim - 1, tdim); - if (!connectivity) - { - throw std::runtime_error( - "Facet-cell connectivity has not been computed."); - } - std::vector cells; - for (auto f : facets) - { - auto c = connectivity->links(f); - assert(c.size() == 1); - cells.push_back(c[0]); - } - sparsitybuild::cells(*sp, cells, dofmaps); - } - */ + 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, int>>, @@ -248,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 Mat create_matrix_nest( const std::vector*>>& a, @@ -265,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 mats(rows * cols, nullptr); + std::shared_ptr> mesh; for (int i = 0; i < rows; ++i) { for (int j = 0; j < cols; ++j) { if (const Form* 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);