DOLFIN-X
DOLFIN-X C++ interface
assembler.h
1 // Copyright (C) 2018-2020 Garth N. Wells
2 //
3 // This file is part of DOLFINX (https://www.fenicsproject.org)
4 //
5 // SPDX-License-Identifier: LGPL-3.0-or-later
6 
7 #pragma once
8 
9 #include "assemble_matrix_impl.h"
10 #include "assemble_scalar_impl.h"
11 #include "assemble_vector_impl.h"
12 #include <Eigen/Dense>
13 #include <Eigen/Sparse>
14 #include <memory>
15 #include <vector>
16 
17 namespace dolfinx
18 {
19 namespace function
20 {
21 class FunctionSpace;
22 } // namespace function
23 
24 namespace fem
25 {
26 template <typename T>
28 template <typename T>
29 class Form;
30 
31 // -- Scalar ----------------------------------------------------------------
32 
38 template <typename T>
40 {
41  return fem::impl::assemble_scalar(M);
42 }
43 
44 // -- Vectors ----------------------------------------------------------------
45 
50 template <typename T>
51 void assemble_vector(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
52  const Form<T>& L)
53 {
54  fem::impl::assemble_vector(b, L);
55 }
56 
57 // FIXME: clarify how x0 is used
58 // FIXME: if bcs entries are set
59 
60 // FIXME: need to pass an array of Vec for x0?
61 // FIXME: clarify zeroing of vector
62 
75 template <typename T>
77  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
78  const std::vector<std::shared_ptr<const Form<T>>>& a,
79  const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
80  const std::vector<Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
81  x0,
82  double scale)
83 {
84  fem::impl::apply_lifting(b, a, bcs1, x0, scale);
85 }
86 
87 // -- Matrices ---------------------------------------------------------------
88 
89 // Experimental
95 template <typename T>
97  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
98  const std::int32_t*, const T*)>& mat_add,
99  const Form<T>& a,
100  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs)
101 {
102 
103  // Index maps for dof ranges
104  auto map0 = a.function_space(0)->dofmap()->index_map;
105  auto map1 = a.function_space(1)->dofmap()->index_map;
106 
107  // Build dof markers
108  std::vector<bool> dof_marker0, dof_marker1;
109  std::int32_t dim0
110  = map0->block_size() * (map0->size_local() + map0->num_ghosts());
111  std::int32_t dim1
112  = map1->block_size() * (map1->size_local() + map1->num_ghosts());
113  for (std::size_t k = 0; k < bcs.size(); ++k)
114  {
115  assert(bcs[k]);
116  assert(bcs[k]->function_space());
117  if (a.function_space(0)->contains(*bcs[k]->function_space()))
118  {
119  dof_marker0.resize(dim0, false);
120  bcs[k]->mark_dofs(dof_marker0);
121  }
122  if (a.function_space(1)->contains(*bcs[k]->function_space()))
123  {
124  dof_marker1.resize(dim1, false);
125  bcs[k]->mark_dofs(dof_marker1);
126  }
127  }
128 
129  // Assemble
130  impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
131 }
132 
143 template <typename T>
145  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
146  const std::int32_t*, const T*)>& mat_add,
147  const Form<T>& a, const std::vector<bool>& dof_marker0,
148  const std::vector<bool>& dof_marker1)
149 
150 {
151  impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
152 }
153 
160 template <typename T>
161 Eigen::SparseMatrix<T, Eigen::RowMajor> assemble_matrix_eigen(
162  const Form<T>& a,
163  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs)
164 {
165  // Lambda function creating Eigen::Triplet array
166  std::vector<Eigen::Triplet<T>> triplets;
167  const auto mat_add
168  = [&triplets](std::int32_t nrow, const std::int32_t* rows,
169  std::int32_t ncol, const std::int32_t* cols, const T* v) {
170  for (int i = 0; i < nrow; ++i)
171  for (int j = 0; j < ncol; ++j)
172  triplets.emplace_back(rows[i], cols[j], v[i * ncol + j]);
173  return 0;
174  };
175 
176  // Assemble
177  assemble_matrix<T>(mat_add, a, bcs);
178 
179  auto map0 = a.function_space(0)->dofmap()->index_map;
180  auto map1 = a.function_space(1)->dofmap()->index_map;
181  Eigen::SparseMatrix<T, Eigen::RowMajor> mat(
182  map0->block_size() * (map0->size_local() + map0->num_ghosts()),
183  map1->block_size() * (map1->size_local() + map1->num_ghosts()));
184  mat.setFromTriplets(triplets.begin(), triplets.end());
185  return mat;
186 }
187 
198 template <typename T>
200  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
201  const std::int32_t*, const T*)>& mat_add,
202  const Eigen::Ref<const Eigen::Array<std::int32_t, Eigen::Dynamic, 1>>& rows,
203  T diagonal = 1.0)
204 {
205  for (Eigen::Index i = 0; i < rows.size(); ++i)
206  {
207  const std::int32_t row = rows(i);
208  mat_add(1, &row, 1, &row, &diagonal);
209  }
210 }
211 
227 template <typename T>
229  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
230  const std::int32_t*, const T*)>& mat_add,
231  const function::FunctionSpace& V,
232  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
233  T diagonal = 1.0)
234 {
235  for (const auto& bc : bcs)
236  {
237  assert(bc);
238  if (V.contains(*bc->function_space()))
239  add_diagonal<T>(mat_add, bc->dofs_owned().col(0), diagonal);
240  }
241 }
242 
243 // -- Setting bcs ------------------------------------------------------------
244 
245 // FIXME: Move these function elsewhere?
246 
247 // FIXME: clarify x0
248 // FIXME: clarify what happens with ghosts
249 
253 template <typename T>
254 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
255  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
256  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
257  double scale = 1.0)
258 {
259  if (b.rows() > x0.rows())
260  throw std::runtime_error("Size mismatch between b and x0 vectors.");
261  for (const auto& bc : bcs)
262  {
263  assert(bc);
264  bc->set(b, x0, scale);
265  }
266 }
267 
271 template <typename T>
272 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
273  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
274  double scale = 1.0)
275 {
276  for (const auto& bc : bcs)
277  {
278  assert(bc);
279  bc->set(b, scale);
280  }
281 }
282 
283 // FIXME: Handle null block
284 // FIXME: Pass function spaces rather than forms
285 
293 template <typename T>
294 std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>
295 bcs_rows(const std::vector<const Form<T>*>& L,
296  const std::vector<std::shared_ptr<const fem::DirichletBC<T>>>& bcs)
297 {
298  // Pack DirichletBC pointers for rows
299  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>> bcs0(
300  L.size());
301  for (std::size_t i = 0; i < L.size(); ++i)
302  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs)
303  if (L[i]->function_space(0)->contains(*bc->function_space()))
304  bcs0[i].push_back(bc);
305  return bcs0;
306 }
307 
308 // FIXME: Handle null block
309 // FIXME: Pass function spaces rather than forms
310 
318 template <typename T>
319 std::vector<
320  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
321 bcs_cols(const std::vector<std::vector<std::shared_ptr<const Form<T>>>>& a,
322  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs)
323 {
324  // Pack DirichletBC pointers for columns
325  std::vector<
326  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
327  bcs1(a.size());
328  for (std::size_t i = 0; i < a.size(); ++i)
329  {
330  for (std::size_t j = 0; j < a[i].size(); ++j)
331  {
332  bcs1[i].resize(a[j].size());
333  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs)
334  {
335  // FIXME: handle case where a[i][j] is null
336  if (a[i][j])
337  {
338  if (a[i][j]->function_space(1)->contains(*bc->function_space()))
339  bcs1[i][j].push_back(bc);
340  }
341  }
342  }
343  }
344 
345  return bcs1;
346 }
347 
348 } // namespace fem
349 } // namespace dolfinx
dolfinx::fem::apply_lifting
void apply_lifting(Eigen::Ref< Eigen::Matrix< T, Eigen::Dynamic, 1 >> b, const std::vector< std::shared_ptr< const Form< T >>> &a, const std::vector< std::vector< std::shared_ptr< const DirichletBC< T >>>> &bcs1, const std::vector< Eigen::Ref< const Eigen::Matrix< T, Eigen::Dynamic, 1 >>> &x0, double scale)
Modify b such that:
Definition: assembler.h:76
dolfinx::fem::assemble_matrix_eigen
Eigen::SparseMatrix< T, Eigen::RowMajor > assemble_matrix_eigen(const Form< T > &a, const std::vector< std::shared_ptr< const DirichletBC< T >>> &bcs)
Definition: assembler.h:161
dolfinx::fem::Form
Class for variational forms.
Definition: assembler.h:29
dolfinx::fem::set_bc
void set_bc(Eigen::Ref< Eigen::Matrix< T, Eigen::Dynamic, 1 >> b, const std::vector< std::shared_ptr< const DirichletBC< T >>> &bcs, const Eigen::Ref< const Eigen::Matrix< T, Eigen::Dynamic, 1 >> &x0, double scale=1.0)
Set bc values in owned (local) part of the PETScVector, multiplied by 'scale'. The vectors b and x0 m...
Definition: assembler.h:254
dolfinx::fem::assemble_vector
void assemble_vector(Eigen::Ref< Eigen::Matrix< T, Eigen::Dynamic, 1 >> b, const Form< T > &L)
Assemble linear form into an Eigen vector.
Definition: assembler.h:51
dolfinx::fem::add_diagonal
void add_diagonal(const std::function< int(std::int32_t, const std::int32_t *, std::int32_t, const std::int32_t *, const T *)> &mat_add, const Eigen::Ref< const Eigen::Array< std::int32_t, Eigen::Dynamic, 1 >> &rows, T diagonal=1.0)
Adds a value to the diagonal of a matrix for specified rows. It is typically called after assembly....
Definition: assembler.h:199
dolfinx::fem::assemble_scalar
T assemble_scalar(const Form< T > &M)
Assemble functional into scalar. Caller is responsible for accumulation across processes.
Definition: assembler.h:39
dolfinx::fem::Form::function_space
std::shared_ptr< const function::FunctionSpace > function_space(int i) const
Return function space for given argument.
Definition: Form.h:235
dolfinx::fem::DirichletBC
Interface for setting (strong) Dirichlet boundary conditions.
Definition: assembler.h:27
dolfinx::fem::bcs_cols
std::vector< std::vector< std::vector< std::shared_ptr< const fem::DirichletBC< T > > > > > bcs_cols(const std::vector< std::vector< std::shared_ptr< const Form< T >>>> &a, const std::vector< std::shared_ptr< const DirichletBC< T >>> &bcs)
Arrange boundary conditions by block.
Definition: assembler.h:321
dolfinx::function::FunctionSpace::contains
bool contains(const FunctionSpace &V) const
Check whether V is subspace of this, or this itself.
Definition: FunctionSpace.cpp:252
dolfinx::function::FunctionSpace
This class represents a finite element function space defined by a mesh, a finite element,...
Definition: FunctionSpace.h:38
dolfinx::fem::bcs_rows
std::vector< std::vector< std::shared_ptr< const fem::DirichletBC< T > > > > bcs_rows(const std::vector< const Form< T > * > &L, const std::vector< std::shared_ptr< const fem::DirichletBC< T >>> &bcs)
Arrange boundary conditions by block.
Definition: assembler.h:295
dolfinx::fem::assemble_matrix
void assemble_matrix(const std::function< int(std::int32_t, const std::int32_t *, std::int32_t, const std::int32_t *, const T *)> &mat_add, const Form< T > &a, const std::vector< std::shared_ptr< const DirichletBC< T >>> &bcs)
Assemble bilinear form into a matrix.
Definition: assembler.h:96