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 <memory>
14 #include <vector>
15 
16 namespace dolfinx
17 {
18 namespace function
19 {
20 class FunctionSpace;
21 } // namespace function
22 
23 namespace fem
24 {
25 template <typename T>
27 template <typename T>
28 class Form;
29 
30 // -- Scalar ----------------------------------------------------------------
31 
37 template <typename T>
39 {
40  return fem::impl::assemble_scalar(M);
41 }
42 
43 // -- Vectors ----------------------------------------------------------------
44 
49 template <typename T>
50 void assemble_vector(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
51  const Form<T>& L)
52 {
53  fem::impl::assemble_vector(b, L);
54 }
55 
56 // FIXME: clarify how x0 is used
57 // FIXME: if bcs entries are set
58 
59 // FIXME: need to pass an array of Vec for x0?
60 // FIXME: clarify zeroing of vector
61 
74 template <typename T>
76  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
77  const std::vector<std::shared_ptr<const Form<T>>>& a,
78  const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
79  const std::vector<Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
80  x0,
81  double scale)
82 {
83  fem::impl::apply_lifting(b, a, bcs1, x0, scale);
84 }
85 
86 // -- Matrices ---------------------------------------------------------------
87 
88 // Experimental
94 template <typename T>
96  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
97  const std::int32_t*, const T*)>& mat_add,
98  const Form<T>& a,
99  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs)
100 {
101 
102  // Index maps for dof ranges
103  auto map0 = a.function_space(0)->dofmap()->index_map;
104  auto map1 = a.function_space(1)->dofmap()->index_map;
105 
106  // Build dof markers
107  std::vector<bool> dof_marker0, dof_marker1;
108  std::int32_t dim0
109  = map0->block_size() * (map0->size_local() + map0->num_ghosts());
110  std::int32_t dim1
111  = map1->block_size() * (map1->size_local() + map1->num_ghosts());
112  for (std::size_t k = 0; k < bcs.size(); ++k)
113  {
114  assert(bcs[k]);
115  assert(bcs[k]->function_space());
116  if (a.function_space(0)->contains(*bcs[k]->function_space()))
117  {
118  dof_marker0.resize(dim0, false);
119  bcs[k]->mark_dofs(dof_marker0);
120  }
121  if (a.function_space(1)->contains(*bcs[k]->function_space()))
122  {
123  dof_marker1.resize(dim1, false);
124  bcs[k]->mark_dofs(dof_marker1);
125  }
126  }
127 
128  // Assemble
129  impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
130 }
131 
142 template <typename T>
144  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
145  const std::int32_t*, const T*)>& mat_add,
146  const Form<T>& a, const std::vector<bool>& dof_marker0,
147  const std::vector<bool>& dof_marker1)
148 
149 {
150  impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
151 }
152 
163 template <typename T>
165  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
166  const std::int32_t*, const T*)>& mat_add,
167  const Eigen::Ref<const Eigen::Array<std::int32_t, Eigen::Dynamic, 1>>& rows,
168  T diagonal = 1.0)
169 {
170  for (Eigen::Index i = 0; i < rows.size(); ++i)
171  {
172  const std::int32_t row = rows(i);
173  mat_add(1, &row, 1, &row, &diagonal);
174  }
175 }
176 
192 template <typename T>
194  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
195  const std::int32_t*, const T*)>& mat_add,
196  const function::FunctionSpace& V,
197  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
198  T diagonal = 1.0)
199 {
200  for (const auto& bc : bcs)
201  {
202  assert(bc);
203  if (V.contains(*bc->function_space()))
204  add_diagonal<T>(mat_add, bc->dofs_owned().col(0), diagonal);
205  }
206 }
207 
208 // -- Setting bcs ------------------------------------------------------------
209 
210 // FIXME: Move these function elsewhere?
211 
212 // FIXME: clarify x0
213 // FIXME: clarify what happens with ghosts
214 
218 template <typename T>
219 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
220  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
221  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
222  double scale = 1.0)
223 {
224  if (b.rows() > x0.rows())
225  throw std::runtime_error("Size mismatch between b and x0 vectors.");
226  for (const auto& bc : bcs)
227  {
228  assert(bc);
229  bc->set(b, x0, scale);
230  }
231 }
232 
236 template <typename T>
237 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
238  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
239  double scale = 1.0)
240 {
241  for (const auto& bc : bcs)
242  {
243  assert(bc);
244  bc->set(b, scale);
245  }
246 }
247 
248 // FIXME: Handle null block
249 // FIXME: Pass function spaces rather than forms
250 
258 template <typename T>
259 std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>
260 bcs_rows(const std::vector<const Form<T>*>& L,
261  const std::vector<std::shared_ptr<const fem::DirichletBC<T>>>& bcs)
262 {
263  // Pack DirichletBC pointers for rows
264  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>> bcs0(
265  L.size());
266  for (std::size_t i = 0; i < L.size(); ++i)
267  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs)
268  if (L[i]->function_space(0)->contains(*bc->function_space()))
269  bcs0[i].push_back(bc);
270  return bcs0;
271 }
272 
273 // FIXME: Handle null block
274 // FIXME: Pass function spaces rather than forms
275 
283 template <typename T>
284 std::vector<
285  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
286 bcs_cols(const std::vector<std::vector<std::shared_ptr<const Form<T>>>>& a,
287  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs)
288 {
289  // Pack DirichletBC pointers for columns
290  std::vector<
291  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
292  bcs1(a.size());
293  for (std::size_t i = 0; i < a.size(); ++i)
294  {
295  for (std::size_t j = 0; j < a[i].size(); ++j)
296  {
297  bcs1[i].resize(a[j].size());
298  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs)
299  {
300  // FIXME: handle case where a[i][j] is null
301  if (a[i][j])
302  {
303  if (a[i][j]->function_space(1)->contains(*bc->function_space()))
304  bcs1[i][j].push_back(bc);
305  }
306  }
307  }
308  }
309 
310  return bcs1;
311 }
312 
313 } // namespace fem
314 } // 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:75
dolfinx::fem::Form
Class for variational forms.
Definition: assembler.h:28
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:219
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:50
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:164
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:38
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:204
dolfinx::fem::DirichletBC
Interface for setting (strong) Dirichlet boundary conditions.
Definition: assembler.h:26
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:286
dolfinx::function::FunctionSpace::contains
bool contains(const FunctionSpace &V) const
Check whether V is subspace of this, or this itself.
Definition: FunctionSpace.cpp:255
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:260
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:95