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>
26 class DirichletBC;
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_spaces().at(0)->dofmap()->index_map;
104  auto map1 = a.function_spaces().at(1)->dofmap()->index_map;
105 
106  // Build dof markers
107  std::vector<bool> dof_marker0, dof_marker1;
108  assert(map0);
109  std::int32_t dim0
110  = map0->block_size() * (map0->size_local() + map0->num_ghosts());
111  assert(map1);
112  std::int32_t dim1
113  = map1->block_size() * (map1->size_local() + map1->num_ghosts());
114  for (std::size_t k = 0; k < bcs.size(); ++k)
115  {
116  assert(bcs[k]);
117  assert(bcs[k]->function_space());
118  if (a.function_spaces().at(0)->contains(*bcs[k]->function_space()))
119  {
120  dof_marker0.resize(dim0, false);
121  bcs[k]->mark_dofs(dof_marker0);
122  }
123  if (a.function_spaces().at(1)->contains(*bcs[k]->function_space()))
124  {
125  dof_marker1.resize(dim1, false);
126  bcs[k]->mark_dofs(dof_marker1);
127  }
128  }
129 
130  // Assemble
131  impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
132 }
133 
144 template <typename T>
146  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
147  const std::int32_t*, const T*)>& mat_add,
148  const Form<T>& a, const std::vector<bool>& dof_marker0,
149  const std::vector<bool>& dof_marker1)
150 
151 {
152  impl::assemble_matrix(mat_add, a, dof_marker0, dof_marker1);
153 }
154 
165 template <typename T>
167  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
168  const std::int32_t*, const T*)>& mat_add,
169  const Eigen::Ref<const Eigen::Array<std::int32_t, Eigen::Dynamic, 1>>& rows,
170  T diagonal = 1.0)
171 {
172  for (Eigen::Index i = 0; i < rows.size(); ++i)
173  {
174  const std::int32_t row = rows(i);
175  mat_add(1, &row, 1, &row, &diagonal);
176  }
177 }
178 
194 template <typename T>
196  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
197  const std::int32_t*, const T*)>& mat_add,
198  const function::FunctionSpace& V,
199  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
200  T diagonal = 1.0)
201 {
202  for (const auto& bc : bcs)
203  {
204  assert(bc);
205  if (V.contains(*bc->function_space()))
206  add_diagonal<T>(mat_add, bc->dofs_owned().col(0), diagonal);
207  }
208 }
209 
210 // -- Setting bcs ------------------------------------------------------------
211 
212 // FIXME: Move these function elsewhere?
213 
214 // FIXME: clarify x0
215 // FIXME: clarify what happens with ghosts
216 
220 template <typename T>
221 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
222  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
223  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
224  double scale = 1.0)
225 {
226  if (b.rows() > x0.rows())
227  throw std::runtime_error("Size mismatch between b and x0 vectors.");
228  for (const auto& bc : bcs)
229  {
230  assert(bc);
231  bc->set(b, x0, scale);
232  }
233 }
234 
238 template <typename T>
239 void set_bc(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
240  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs,
241  double scale = 1.0)
242 {
243  for (const auto& bc : bcs)
244  {
245  assert(bc);
246  bc->set(b, scale);
247  }
248 }
249 
250 // FIXME: Handle null block
251 // FIXME: Pass function spaces rather than forms
252 
260 template <typename T>
261 std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>
262 bcs_rows(const std::vector<const Form<T>*>& L,
263  const std::vector<std::shared_ptr<const fem::DirichletBC<T>>>& bcs)
264 {
265  // Pack DirichletBC pointers for rows
266  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>> bcs0(
267  L.size());
268  for (std::size_t i = 0; i < L.size(); ++i)
269  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs)
270  if (L[i]->function_spaces()[0]->contains(*bc->function_space()))
271  bcs0[i].push_back(bc);
272  return bcs0;
273 }
274 
275 // FIXME: Handle null block
276 // FIXME: Pass function spaces rather than forms
277 
285 template <typename T>
286 std::vector<
287  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
288 bcs_cols(const std::vector<std::vector<std::shared_ptr<const Form<T>>>>& a,
289  const std::vector<std::shared_ptr<const DirichletBC<T>>>& bcs)
290 {
291  // Pack DirichletBC pointers for columns
292  std::vector<
293  std::vector<std::vector<std::shared_ptr<const fem::DirichletBC<T>>>>>
294  bcs1(a.size());
295  for (std::size_t i = 0; i < a.size(); ++i)
296  {
297  for (std::size_t j = 0; j < a[i].size(); ++j)
298  {
299  bcs1[i].resize(a[j].size());
300  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs)
301  {
302  // FIXME: handle case where a[i][j] is null
303  if (a[i][j])
304  {
305  if (a[i][j]->function_spaces()[1]->contains(*bc->function_space()))
306  bcs1[i][j].push_back(bc);
307  }
308  }
309  }
310  }
311 
312  return bcs1;
313 }
314 
315 } // namespace fem
316 } // 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::function_spaces
const std::vector< std::shared_ptr< const function::FunctionSpace > > & function_spaces() const
Return function spaces for all arguments.
Definition: Form.h:157
dolfinx::fem::Form
Class for variational forms.
Definition: Form.h:66
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:221
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:166
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::DirichletBC
Interface for setting (strong) Dirichlet boundary conditions.
Definition: DirichletBC.h:102
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:288
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:262
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