13#include <dolfinx/la/petsc.h>
29template <dolfinx::scalar T, std::
floating_po
int U>
41template <std::
floating_po
int T>
43 std::string type = std::string())
58template <std::
floating_po
int T>
61 std::string type = std::string())
64 std::array<std::vector<std::shared_ptr<const FunctionSpace<T>>>, 2> V
66 std::array<std::vector<int>, 2> bs_dofs;
67 for (std::size_t i = 0; i < 2; ++i)
70 bs_dofs[i].push_back(_V->dofmap()->bs());
74 std::shared_ptr<const mesh::Mesh<T>> mesh;
75 std::vector<std::vector<std::unique_ptr<la::SparsityPattern>>> patterns(
77 for (std::size_t row = 0; row < V[0].size(); ++row)
79 for (std::size_t col = 0; col < V[1].size(); ++col)
83 patterns[row].push_back(std::make_unique<la::SparsityPattern>(
89 patterns[row].push_back(
nullptr);
94 throw std::runtime_error(
"Could not find a Mesh.");
97 std::array<std::vector<std::pair<
98 std::reference_wrapper<const common::IndexMap>,
int>>,
101 for (std::size_t d = 0; d < 2; ++d)
103 for (
auto& space : V[d])
105 maps[d].emplace_back(*space->dofmap()->index_map,
106 space->dofmap()->index_map_bs());
111 std::vector<std::vector<const la::SparsityPattern*>> p(V[0].size());
112 for (std::size_t row = 0; row < V[0].size(); ++row)
113 for (std::size_t col = 0; col < V[1].size(); ++col)
114 p[row].push_back(patterns[row][col].get());
127 std::array<std::vector<PetscInt>, 2> _maps;
128 for (
int d = 0; d < 2; ++d)
138 std::pair<std::reference_wrapper<const common::IndexMap>,
int>>& map
140 std::vector<PetscInt>& _map = _maps[d];
143 const auto [rank_offset, local_offset, ghosts, _]
145 for (std::size_t f = 0; f < map.size(); ++f)
147 auto offset = local_offset[f];
149 int bs = map[f].second;
150 for (std::int32_t i = 0; i < bs * imap.
size_local(); ++i)
151 _map.push_back(i + rank_offset + offset);
152 for (std::int32_t i = 0; i < bs * imap.
num_ghosts(); ++i)
153 _map.push_back(ghosts[f][i]);
158 ISLocalToGlobalMapping petsc_local_to_global0;
159 ISLocalToGlobalMappingCreate(MPI_COMM_SELF, 1, _maps[0].size(),
160 _maps[0].data(), PETSC_COPY_VALUES,
161 &petsc_local_to_global0);
164 MatSetLocalToGlobalMapping(A, petsc_local_to_global0,
165 petsc_local_to_global0);
166 ISLocalToGlobalMappingDestroy(&petsc_local_to_global0);
171 ISLocalToGlobalMapping petsc_local_to_global1;
172 ISLocalToGlobalMappingCreate(MPI_COMM_SELF, 1, _maps[1].size(),
173 _maps[1].data(), PETSC_COPY_VALUES,
174 &petsc_local_to_global1);
175 MatSetLocalToGlobalMapping(A, petsc_local_to_global0,
176 petsc_local_to_global1);
177 ISLocalToGlobalMappingDestroy(&petsc_local_to_global0);
178 ISLocalToGlobalMappingDestroy(&petsc_local_to_global1);
187template <std::
floating_po
int T>
190 const std::vector<std::vector<std::string>>& types)
195 std::vector<std::vector<std::string>> _types(
196 a.size(), std::vector<std::string>(a.front().size()));
202 int cols = a.front().size();
203 std::vector<Mat> mats(rows * cols,
nullptr);
204 std::shared_ptr<const mesh::Mesh<T>> mesh;
205 for (
int i = 0; i < rows; ++i)
207 for (
int j = 0; j < cols; ++j)
218 throw std::runtime_error(
"Could not find a Mesh.");
222 MatCreate(mesh->comm(), &A);
223 MatSetType(A, MATNEST);
224 MatNestSetSubMats(A, rows,
nullptr, cols,
nullptr, mats.data());
228 for (std::size_t i = 0; i < mats.size(); ++i)
231 MatDestroy(&mats[i]);
242 std::pair<std::reference_wrapper<const common::IndexMap>,
int>>& maps);
247 std::pair<std::reference_wrapper<const common::IndexMap>,
int>>& maps);
263template <std::
floating_po
int T>
266 std::span<const PetscScalar> constants,
267 const std::map<std::pair<IntegralType, int>,
268 std::pair<std::span<const PetscScalar>,
int>>& coeffs)
271 VecGhostGetLocalForm(b, &b_local);
273 VecGetSize(b_local, &n);
274 PetscScalar* array =
nullptr;
275 VecGetArray(b_local, &array);
276 std::span<PetscScalar> _b(array, n);
278 VecRestoreArray(b_local, &array);
279 VecGhostRestoreLocalForm(b, &b_local);
291template <std::
floating_po
int T>
295 VecGhostGetLocalForm(b, &b_local);
297 VecGetSize(b_local, &n);
298 PetscScalar* array =
nullptr;
299 VecGetArray(b_local, &array);
300 std::span<PetscScalar> _b(array, n);
302 VecRestoreArray(b_local, &array);
303 VecGhostRestoreLocalForm(b, &b_local);
327template <std::
floating_po
int T>
330 const std::vector<std::span<const PetscScalar>>& constants,
331 const std::vector<std::map<std::pair<IntegralType, int>,
332 std::pair<std::span<const PetscScalar>,
int>>>&
336 const std::vector<Vec>& x0, PetscScalar scale)
339 VecGhostGetLocalForm(b, &b_local);
341 VecGetSize(b_local, &n);
342 PetscScalar* array =
nullptr;
343 VecGetArray(b_local, &array);
344 std::span<PetscScalar> _b(array, n);
350 std::vector<std::span<const PetscScalar>> x0_ref;
351 std::vector<Vec> x0_local(a.size());
352 std::vector<const PetscScalar*> x0_array(a.size());
353 for (std::size_t i = 0; i < a.size(); ++i)
356 VecGhostGetLocalForm(x0[i], &x0_local[i]);
358 VecGetSize(x0_local[i], &n);
359 VecGetArrayRead(x0_local[i], &x0_array[i]);
360 x0_ref.emplace_back(x0_array[i], n);
363 std::vector x0_tmp(x0_ref.begin(), x0_ref.end());
366 for (std::size_t i = 0; i < x0_local.size(); ++i)
368 VecRestoreArrayRead(x0_local[i], &x0_array[i]);
369 VecGhostRestoreLocalForm(x0[i], &x0_local[i]);
373 VecRestoreArray(b_local, &array);
374 VecGhostRestoreLocalForm(b, &b_local);
395template <std::
floating_po
int T>
402 const std::vector<Vec>& x0, PetscScalar scale)
405 VecGhostGetLocalForm(b, &b_local);
407 VecGetSize(b_local, &n);
408 PetscScalar* array =
nullptr;
409 VecGetArray(b_local, &array);
410 std::span<PetscScalar> _b(array, n);
413 fem::apply_lifting<PetscScalar>(_b, a, bcs1, {}, scale);
416 std::vector<std::span<const PetscScalar>> x0_ref;
417 std::vector<Vec> x0_local(a.size());
418 std::vector<const PetscScalar*> x0_array(a.size());
419 for (std::size_t i = 0; i < a.size(); ++i)
422 VecGhostGetLocalForm(x0[i], &x0_local[i]);
424 VecGetSize(x0_local[i], &n);
425 VecGetArrayRead(x0_local[i], &x0_array[i]);
426 x0_ref.emplace_back(x0_array[i], n);
429 std::vector x0_tmp(x0_ref.begin(), x0_ref.end());
430 fem::apply_lifting<PetscScalar>(_b, a, bcs1, x0_tmp, scale);
432 for (std::size_t i = 0; i < x0_local.size(); ++i)
434 VecRestoreArrayRead(x0_local[i], &x0_array[i]);
435 VecGhostRestoreLocalForm(x0[i], &x0_local[i]);
439 VecRestoreArray(b_local, &array);
440 VecGhostRestoreLocalForm(b, &b_local);
453template <std::
floating_po
int T>
457 const Vec x0, PetscScalar scale = 1)
460 VecGetLocalSize(b, &n);
461 PetscScalar* array =
nullptr;
462 VecGetArray(b, &array);
463 std::span<PetscScalar> _b(array, n);
467 VecGhostGetLocalForm(x0, &x0_local);
469 VecGetSize(x0_local, &n);
470 const PetscScalar* array =
nullptr;
471 VecGetArrayRead(x0_local, &array);
472 std::span<const PetscScalar> _x0(array, n);
474 VecRestoreArrayRead(x0_local, &array);
475 VecGhostRestoreLocalForm(x0, &x0_local);
480 VecRestoreArray(b, &array);
This class represents the distribution index arrays across processes. An index array is a contiguous ...
Definition IndexMap.h:94
std::int32_t num_ghosts() const noexcept
Number of ghost indices on this process.
Definition IndexMap.cpp:866
std::int32_t size_local() const noexcept
Number of indices owned by this process.
Definition IndexMap.cpp:868
Object for setting (strong) Dirichlet boundary conditions.
Definition DirichletBC.h:259
Sparsity pattern data structure that can be used to initialize sparse matrices. After assembly,...
Definition SparsityPattern.h:26
void finalize()
Finalize sparsity pattern and communicate off-process entries.
Definition SparsityPattern.cpp:226
Miscellaneous classes, functions and types.
Definition dolfinx_common.h:8
std::tuple< std::int64_t, std::vector< std::int32_t >, std::vector< std::vector< std::int64_t > >, std::vector< std::vector< int > > > stack_index_maps(const std::vector< std::pair< std::reference_wrapper< const IndexMap >, int > > &maps)
Compute layout data and ghost indices for a stacked (concatenated) index map, i.e....
Definition IndexMap.cpp:583
void set_bc(Vec b, const std::vector< std::shared_ptr< const DirichletBC< PetscScalar, T > > > &bcs, const Vec x0, PetscScalar scale=1)
Set bc values in owned (local) part of the PETSc vector, multiplied by 'scale'. The vectors b and x0 ...
Definition petsc.h:454
void apply_lifting(Vec b, const std::vector< std::shared_ptr< const Form< PetscScalar, T > > > &a, const std::vector< std::span< const PetscScalar > > &constants, const std::vector< std::map< std::pair< IntegralType, int >, std::pair< std::span< const PetscScalar >, int > > > &coeffs, const std::vector< std::vector< std::shared_ptr< const DirichletBC< PetscScalar, T > > > > &bcs1, const std::vector< Vec > &x0, PetscScalar scale)
Modify RHS vector to account for Dirichlet boundary conditions.
Definition petsc.h:328
Mat create_matrix(const Form< PetscScalar, T > &a, std::string type=std::string())
Create a matrix.
Definition petsc.h:42
void assemble_vector(Vec b, const Form< PetscScalar, T > &L, std::span< const PetscScalar > constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const PetscScalar >, int > > &coeffs)
Assemble linear form into an already allocated PETSc vector.
Definition petsc.h:264
Mat create_matrix_block(const std::vector< std::vector< const Form< PetscScalar, T > * > > &a, std::string type=std::string())
Initialise a monolithic matrix for an array of bilinear forms.
Definition petsc.h:59
Mat create_matrix_nest(const std::vector< std::vector< const Form< PetscScalar, T > * > > &a, const std::vector< std::vector< std::string > > &types)
Create nested (MatNest) matrix.
Definition petsc.h:188
Vec create_vector_block(const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Initialise monolithic vector. Vector is not zeroed.
Definition petsc.cpp:15
Vec create_vector_nest(const std::vector< std::pair< std::reference_wrapper< const common::IndexMap >, int > > &maps)
Create nested (VecNest) vector. Vector is not zeroed.
Definition petsc.cpp:60
Finite element method functionality.
Definition assemble_matrix_impl.h:26
void assemble_vector(std::span< T > b, const Form< T, U > &L, std::span< const T > constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients)
Assemble linear form into a vector.
Definition assembler.h:109
void apply_lifting(std::span< T > b, const std::vector< std::shared_ptr< const Form< T, U > > > &a, const std::vector< std::span< const T > > &constants, const std::vector< std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > > &coeffs, const std::vector< std::vector< std::shared_ptr< const DirichletBC< T, U > > > > &bcs1, const std::vector< std::span< const T > > &x0, T scale)
Modify b such that:
Definition assembler.h:150
std::vector< std::vector< std::array< std::shared_ptr< const FunctionSpace< U > >, 2 > > > extract_function_spaces(const std::vector< std::vector< const Form< T, U > * > > &a)
Extract test (0) and trial (1) function spaces pairs for each bilinear form for a rectangular array o...
Definition utils.h:125
void set_bc(std::span< T > b, const std::vector< std::shared_ptr< const DirichletBC< T, U > > > &bcs, std::span< const T > x0, T scale=1)
Set bc values in owned (local) part of the vector, multiplied by 'scale'. The vectors b and x0 must h...
Definition assembler.h:438
la::SparsityPattern create_sparsity_pattern(const Form< T, U > &a)
Create a sparsity pattern for a given form.
Definition utils.h:150
std::array< std::vector< std::shared_ptr< const FunctionSpace< T > > >, 2 > common_function_spaces(const std::vector< std::vector< std::array< std::shared_ptr< const FunctionSpace< T > >, 2 > > > &V)
Extract FunctionSpaces for (0) rows blocks and (1) columns blocks from a rectangular array of (test,...
Definition FunctionSpace.h:374
Mat create_matrix(MPI_Comm comm, const SparsityPattern &sp, std::string type=std::string())
Create a PETSc Mat. Caller is responsible for destroying the returned object.
Definition petsc.cpp:232