DOLFIN-X
DOLFIN-X C++ interface
assemble_vector_impl.h
1 // Copyright (C) 2018-2019 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 "DirichletBC.h"
10 #include "DofMap.h"
11 #include "Form.h"
12 #include "utils.h"
13 #include <dolfinx/common/IndexMap.h>
14 #include <dolfinx/common/types.h>
15 #include <dolfinx/fem/Constant.h>
16 #include <dolfinx/fem/FunctionSpace.h>
17 #include <dolfinx/graph/AdjacencyList.h>
18 #include <dolfinx/mesh/Geometry.h>
19 #include <dolfinx/mesh/Mesh.h>
20 #include <dolfinx/mesh/Topology.h>
21 #include <functional>
22 #include <memory>
23 #include <vector>
24 
25 namespace dolfinx::fem::impl
26 {
27 
29 
34 template <typename T>
35 void assemble_vector(tcb::span<T> b, const Form<T>& L);
36 
38 template <typename T>
39 void assemble_cells(
40  tcb::span<T> b, const mesh::Geometry& geometry,
41  const std::vector<std::int32_t>& active_cells,
42  const graph::AdjacencyList<std::int32_t>& dofmap, const int bs,
43  const std::function<void(T*, const T*, const T*, const double*, const int*,
44  const std::uint8_t*, const std::uint32_t)>& kernel,
45  const array2d<T>& coeffs, const std::vector<T>& constant_values,
46  const std::vector<std::uint32_t>& cell_info);
47 
49 template <typename T>
50 void assemble_exterior_facets(
51  tcb::span<T> b, const mesh::Mesh& mesh,
52  const std::vector<std::int32_t>& active_facets,
53  const graph::AdjacencyList<std::int32_t>& dofmap, const int bs,
54  const std::function<void(T*, const T*, const T*, const double*, const int*,
55  const std::uint8_t*, const std::uint32_t)>& fn,
56  const array2d<T>& coeffs, const std::vector<T>& constant_values,
57  const std::vector<std::uint32_t>& cell_info,
58  const std::vector<std::uint8_t>& perms);
59 
61 template <typename T>
62 void assemble_interior_facets(
63  tcb::span<T> b, const mesh::Mesh& mesh,
64  const std::vector<std::int32_t>& active_facets, const fem::DofMap& dofmap,
65  const std::function<void(T*, const T*, const T*, const double*, const int*,
66  const std::uint8_t*, const std::uint32_t)>& fn,
67  const array2d<T>& coeffs, const std::vector<int>& offsets,
68  const std::vector<T>& constant_values,
69  const std::vector<std::uint32_t>& cell_info,
70  const std::vector<std::uint8_t>& perms);
71 
89 template <typename T>
90 void apply_lifting(
91  tcb::span<T> b, const std::vector<std::shared_ptr<const Form<T>>> a,
92  const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
93  const tcb::span<const T>& x0, double scale);
94 
107 template <typename T>
108 void lift_bc(tcb::span<T> b, const Form<T>& a,
109  const tcb::span<const T>& bc_values1,
110  const std::vector<bool>& bc_markers1, const tcb::span<const T>& x0,
111  double scale);
112 
113 // Implementation of bc application
114 template <typename T>
115 void _lift_bc_cells(
116  tcb::span<T> b, const mesh::Geometry& geometry,
117  const std::function<void(T*, const T*, const T*, const double*, const int*,
118  const std::uint8_t*, const std::uint32_t)>& kernel,
119  const std::vector<std::int32_t>& active_cells,
120  const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,
121  const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
122  const array2d<T>& coeffs, const std::vector<T>& constant_values,
123  const std::vector<std::uint32_t>& cell_info,
124  const tcb::span<const T>& bc_values1, const std::vector<bool>& bc_markers1,
125  const tcb::span<const T>& x0, double scale)
126 {
127  // Prepare cell geometry
128  const int gdim = geometry.dim();
129  const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
130 
131  // FIXME: Add proper interface for num coordinate dofs
132  const int num_dofs_g = x_dofmap.num_links(0);
133  const array2d<double>& x_g = geometry.x();
134 
135  // Data structures used in bc application
136  std::vector<double> coordinate_dofs(num_dofs_g * gdim);
137  std::vector<T> Ae, be;
138 
139  for (std::int32_t c : active_cells)
140  {
141  // Get dof maps for cell
142  auto dmap1 = dofmap1.links(c);
143 
144  // Check if bc is applied to cell
145  bool has_bc = false;
146  for (std::size_t j = 0; j < dmap1.size(); ++j)
147  {
148  for (int k = 0; k < bs1; ++k)
149  {
150  assert(bs1 * dmap1[j] + k < (int)bc_markers1.size());
151  if (bc_markers1[bs1 * dmap1[j] + k])
152  {
153  has_bc = true;
154  break;
155  }
156  }
157  }
158 
159  if (!has_bc)
160  continue;
161 
162  // Get cell coordinates/geometry
163  auto x_dofs = x_dofmap.links(c);
164  for (std::size_t i = 0; i < x_dofs.size(); ++i)
165  {
166  std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
167  std::next(coordinate_dofs.begin(), i * gdim));
168  }
169 
170  // Size data structure for assembly
171  auto dmap0 = dofmap0.links(c);
172 
173  const int num_rows = bs0 * dmap0.size();
174  const int num_cols = bs1 * dmap1.size();
175 
176  auto coeff_array = coeffs.row(c);
177  Ae.resize(num_rows * num_cols);
178  std::fill(Ae.begin(), Ae.end(), 0);
179  kernel(Ae.data(), coeff_array.data(), constant_values.data(),
180  coordinate_dofs.data(), nullptr, nullptr, cell_info[c]);
181 
182  // Size data structure for assembly
183  be.resize(num_rows);
184  std::fill(be.begin(), be.end(), 0);
185  for (std::size_t j = 0; j < dmap1.size(); ++j)
186  {
187  for (int k = 0; k < bs1; ++k)
188  {
189  const std::int32_t jj = bs1 * dmap1[j] + k;
190  assert(jj < (int)bc_markers1.size());
191  if (bc_markers1[jj])
192  {
193  const T bc = bc_values1[jj];
194  const T _x0 = x0.empty() ? 0.0 : x0[jj];
195  // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
196  for (int m = 0; m < num_rows; ++m)
197  be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
198  }
199  }
200  }
201 
202  for (std::size_t i = 0; i < dmap0.size(); ++i)
203  for (int k = 0; k < bs0; ++k)
204  b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
205  }
206 }
207 //----------------------------------------------------------------------------
208 template <typename T>
209 void _lift_bc_exterior_facets(
210  tcb::span<T> b, const mesh::Mesh& mesh,
211  const std::function<void(T*, const T*, const T*, const double*, const int*,
212  const std::uint8_t*, const std::uint32_t)>& kernel,
213  const std::vector<std::int32_t>& active_facets,
214  const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,
215  const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
216  const array2d<T>& coeffs, const std::vector<T>& constant_values,
217  const std::vector<std::uint32_t>& cell_info,
218  const std::vector<std::uint8_t>& perms,
219  const tcb::span<const T>& bc_values1, const std::vector<bool>& bc_markers1,
220  const tcb::span<const T>& x0, double scale)
221 {
222  const int gdim = mesh.geometry().dim();
223  const int tdim = mesh.topology().dim();
224 
225  // Prepare cell geometry
226  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
227 
228  // FIXME: Add proper interface for num coordinate dofs
229  const int num_dofs_g = x_dofmap.num_links(0);
230  const array2d<double>& x_g = mesh.geometry().x();
231 
232  // Data structures used in bc application
233  std::vector<double> coordinate_dofs(num_dofs_g * gdim);
234  std::vector<T> Ae, be;
235 
236  // Iterate over owned facets
237  const mesh::Topology& topology = mesh.topology();
238  auto connectivity = topology.connectivity(tdim - 1, tdim);
239  assert(connectivity);
240  auto c_to_f = topology.connectivity(tdim, tdim - 1);
241  assert(c_to_f);
242  auto map = topology.index_map(tdim - 1);
243  assert(map);
244 
245  for (std::int32_t f : active_facets)
246  {
247  // Create attached cell
248  assert(connectivity->num_links(f) == 1);
249  const std::int32_t cell = connectivity->links(f)[0];
250 
251  // Get local index of facet with respect to the cell
252  auto facets = c_to_f->links(cell);
253  auto it = std::find(facets.begin(), facets.end(), f);
254  assert(it != facets.end());
255  const int local_facet = std::distance(facets.begin(), it);
256 
257  // Get dof maps for cell
258  auto dmap1 = dofmap1.links(cell);
259 
260  // Check if bc is applied to cell
261  bool has_bc = false;
262  for (std::size_t j = 0; j < dmap1.size(); ++j)
263  {
264  for (int k = 0; k < bs1; ++k)
265  {
266  if (bc_markers1[bs1 * dmap1[j] + k])
267  {
268  has_bc = true;
269  break;
270  }
271  }
272  }
273 
274  if (!has_bc)
275  continue;
276 
277  // Get cell coordinates/geometry
278  auto x_dofs = x_dofmap.links(cell);
279  for (std::size_t i = 0; i < x_dofs.size(); ++i)
280  {
281  std::copy_n(x_g.row(x_dofs[i]).data(), gdim,
282  std::next(coordinate_dofs.begin(), i * gdim));
283  }
284 
285  // Size data structure for assembly
286  auto dmap0 = dofmap0.links(cell);
287 
288  const int num_rows = bs0 * dmap0.size();
289  const int num_cols = bs1 * dmap1.size();
290 
291  auto coeff_array = coeffs.row(cell);
292  Ae.resize(num_rows * num_cols);
293  std::fill(Ae.begin(), Ae.end(), 0);
294  kernel(Ae.data(), coeff_array.data(), constant_values.data(),
295  coordinate_dofs.data(), &local_facet,
296  &perms[cell * facets.size() + local_facet], cell_info[cell]);
297 
298  // Size data structure for assembly
299  be.resize(num_rows);
300  std::fill(be.begin(), be.end(), 0);
301  for (std::size_t j = 0; j < dmap1.size(); ++j)
302  {
303  for (int k = 0; k < bs1; ++k)
304  {
305  const std::int32_t jj = bs1 * dmap1[j] + k;
306  if (bc_markers1[jj])
307  {
308 
309  const T bc = bc_values1[jj];
310  const T _x0 = x0.empty() ? 0.0 : x0[jj];
311  // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
312  for (int m = 0; m < num_rows; ++m)
313  be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
314  }
315  }
316  }
317 
318  for (std::size_t i = 0; i < dmap0.size(); ++i)
319  for (int k = 0; k < bs0; ++k)
320  b[bs0 * dmap0[i] + k] += be[bs0 * i + k];
321  }
322 }
323 //----------------------------------------------------------------------------
324 template <typename T>
325 void _lift_bc_interior_facets(
326  tcb::span<T> b, const mesh::Mesh& mesh,
327  const std::function<void(T*, const T*, const T*, const double*, const int*,
328  const std::uint8_t*, const std::uint32_t)>& kernel,
329  const std::vector<std::int32_t>& active_facets,
330  const graph::AdjacencyList<std::int32_t>& dofmap0, int bs0,
331  const graph::AdjacencyList<std::int32_t>& dofmap1, int bs1,
332  const array2d<T>& coeffs, const std::vector<int>& offsets,
333  const std::vector<T>& constant_values,
334  const std::vector<std::uint32_t>& cell_info,
335  const std::vector<std::uint8_t>& perms,
336  const tcb::span<const T>& bc_values1, const std::vector<bool>& bc_markers1,
337  const tcb::span<const T>& x0, double scale)
338 {
339  const int gdim = mesh.geometry().dim();
340  const int tdim = mesh.topology().dim();
341 
342  // Prepare cell geometry
343  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
344 
345  // FIXME: Add proper interface for num coordinate dofs
346  const int num_dofs_g = x_dofmap.num_links(0);
347  const array2d<double>& x_g = mesh.geometry().x();
348 
349  // Data structures used in assembly
350  std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
351  std::vector<T> coeff_array(2 * offsets.back());
352  assert(offsets.back() == int(coeffs.shape[1]));
353  std::vector<T> Ae, be;
354 
355  // Temporaries for joint dofmaps
356  std::vector<std::int32_t> dmapjoint0, dmapjoint1;
357 
358  const mesh::Topology& topology = mesh.topology();
359  auto connectivity = topology.connectivity(tdim - 1, tdim);
360  assert(connectivity);
361  auto c_to_f = topology.connectivity(tdim, tdim - 1);
362  assert(c_to_f);
363  auto f_to_c = topology.connectivity(tdim - 1, tdim);
364  assert(f_to_c);
365  auto map = topology.index_map(tdim - 1);
366  assert(map);
367 
368  const int offset_g = gdim * num_dofs_g;
369  for (std::int32_t f : active_facets)
370  {
371  // Create attached cells
372  auto cells = f_to_c->links(f);
373  assert(cells.size() == 2);
374 
375  // Get local index of facet with respect to the cell
376  auto facets0 = c_to_f->links(cells[0]);
377  auto it0 = std::find(facets0.begin(), facets0.end(), f);
378  assert(it0 != facets0.end());
379  const int local_facet0 = std::distance(facets0.begin(), it0);
380  auto facets1 = c_to_f->links(cells[1]);
381  auto it1 = std::find(facets1.begin(), facets1.end(), f);
382  assert(it1 != facets1.end());
383  const int local_facet1 = std::distance(facets1.begin(), it1);
384 
385  const std::array local_facet{local_facet0, local_facet1};
386 
387  // Get cell geometry
388  auto x_dofs0 = x_dofmap.links(cells[0]);
389  auto x_dofs1 = x_dofmap.links(cells[1]);
390  for (int i = 0; i < num_dofs_g; ++i)
391  {
392  for (int j = 0; j < gdim; ++j)
393  {
394  coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
395  coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
396  }
397  }
398 
399  // Get dof maps for cells and pack
400  const tcb::span<const std::int32_t> dmap0_cell0 = dofmap0.links(cells[0]);
401  const tcb::span<const std::int32_t> dmap0_cell1 = dofmap0.links(cells[1]);
402  dmapjoint0.resize(dmap0_cell0.size() + dmap0_cell1.size());
403  std::copy(dmap0_cell0.begin(), dmap0_cell0.end(), dmapjoint0.begin());
404  std::copy(dmap0_cell1.begin(), dmap0_cell1.end(),
405  std::next(dmapjoint0.begin(), dmap0_cell0.size()));
406 
407  const tcb::span<const std::int32_t> dmap1_cell0 = dofmap1.links(cells[0]);
408  const tcb::span<const std::int32_t> dmap1_cell1 = dofmap1.links(cells[1]);
409  dmapjoint1.resize(dmap1_cell0.size() + dmap1_cell1.size());
410  std::copy(dmap1_cell0.begin(), dmap1_cell0.end(), dmapjoint1.begin());
411  std::copy(dmap1_cell1.begin(), dmap1_cell1.end(),
412  std::next(dmapjoint1.begin(), dmap1_cell0.size()));
413 
414  // Check if bc is applied to cell0
415  bool has_bc = false;
416  for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
417  {
418  for (int k = 0; k < bs1; ++k)
419  {
420  if (bc_markers1[bs1 * dmap1_cell0[j] + k])
421  {
422  has_bc = true;
423  break;
424  }
425  }
426  }
427 
428  // Check if bc is applied to cell1
429  for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
430  {
431  for (int k = 0; k < bs1; ++k)
432  {
433  if (bc_markers1[bs1 * dmap1_cell1[j] + k])
434  {
435  has_bc = true;
436  break;
437  }
438  }
439  }
440 
441  if (!has_bc)
442  continue;
443 
444  // Layout for the restricted coefficients is flattened
445  // w[coefficient][restriction][dof]
446  const auto coeff_cell0 = coeffs.row(cells[0]);
447  const auto coeff_cell1 = coeffs.row(cells[1]);
448 
449  // Loop over coefficients
450  for (std::size_t i = 0; i < offsets.size() - 1; ++i)
451  {
452  // Loop over entries for coefficient i
453  const int num_entries = offsets[i + 1] - offsets[i];
454  std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
455  std::next(coeff_array.begin(), 2 * offsets[i]));
456  std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
457  std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
458  }
459 
460  const int num_rows = bs0 * dmapjoint0.size();
461  const int num_cols = bs1 * dmapjoint1.size();
462 
463  // Tabulate tensor
464  Ae.resize(num_rows * num_cols);
465  std::fill(Ae.begin(), Ae.end(), 0);
466  const int facets_per_cell = facets0.size();
467  const std::array perm{perms[cells[0] * facets_per_cell + local_facet[0]],
468  perms[cells[1] * facets_per_cell + local_facet[1]]};
469  kernel(Ae.data(), coeff_array.data(), constant_values.data(),
470  coordinate_dofs.data(), local_facet.data(), perm.data(),
471  cell_info[cells[0]]);
472 
473  be.resize(num_rows);
474  std::fill(be.begin(), be.end(), 0);
475 
476  // Compute b = b - A*b for cell0
477  for (std::size_t j = 0; j < dmap1_cell0.size(); ++j)
478  {
479  for (int k = 0; k < bs1; ++k)
480  {
481  const std::int32_t jj = bs1 * dmap1_cell0[j] + k;
482  if (bc_markers1[jj])
483  {
484  const T bc = bc_values1[jj];
485  const T _x0 = x0.empty() ? 0.0 : x0[jj];
486  // be -= Ae.col(bs1 * j + k) * scale * (bc - _x0);
487  for (int m = 0; m < num_rows; ++m)
488  be[m] -= Ae[m * num_cols + bs1 * j + k] * scale * (bc - _x0);
489  }
490  }
491  }
492 
493  // Compute b = b - A*b for cell1
494  const int offset = bs1 * dmap1_cell0.size();
495  for (std::size_t j = 0; j < dmap1_cell1.size(); ++j)
496  {
497  for (int k = 0; k < bs1; ++k)
498  {
499  const std::int32_t jj = bs1 * dmap1_cell1[j] + k;
500  if (bc_markers1[jj])
501  {
502  const T bc = bc_values1[jj];
503  const T _x0 = x0.empty() ? 0.0 : x0[jj];
504  // be -= Ae.col(offset + bs1 * j + k) * scale * (bc - x0[jj]);
505  for (int m = 0; m < num_rows; ++m)
506  {
507  be[m]
508  -= Ae[m * num_cols + offset + bs1 * j + k] * scale * (bc - _x0);
509  }
510  }
511  }
512  }
513 
514  for (std::size_t i = 0; i < dmap0_cell0.size(); ++i)
515  for (int k = 0; k < bs0; ++k)
516  b[bs0 * dmap0_cell0[i] + k] += be[bs0 * i + k];
517 
518  for (std::size_t i = 0; i < dmap0_cell1.size(); ++i)
519  for (int k = 0; k < bs0; ++k)
520  b[bs0 * dmap0_cell1[i] + k] += be[offset + bs0 * i + k];
521  }
522 }
523 //-----------------------------------------------------------------------------
524 template <typename T>
525 void assemble_vector(tcb::span<T> b, const Form<T>& L)
526 {
527  std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
528  assert(mesh);
529  const int tdim = mesh->topology().dim();
530  const std::int32_t num_cells
531  = mesh->topology().connectivity(tdim, 0)->num_nodes();
532 
533  // Get dofmap data
534  assert(L.function_spaces().at(0));
535  std::shared_ptr<const fem::DofMap> dofmap
536  = L.function_spaces().at(0)->dofmap();
537  assert(dofmap);
538  const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
539  const int bs = dofmap->bs();
540 
541  // Prepare constants
542  const std::vector<T> constant_values = pack_constants(L);
543 
544  // Prepare coefficients
545  const array2d<T> coeffs = pack_coefficients(L);
546 
547  const bool needs_permutation_data = L.needs_permutation_data();
548  if (needs_permutation_data)
549  mesh->topology_mutable().create_entity_permutations();
550  const std::vector<std::uint32_t>& cell_info
551  = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
552  : std::vector<std::uint32_t>(num_cells);
553 
554  for (int i : L.integral_ids(IntegralType::cell))
555  {
556  const auto& fn = L.kernel(IntegralType::cell, i);
557  const std::vector<std::int32_t>& active_cells
558  = L.domains(IntegralType::cell, i);
559  fem::impl::assemble_cells(b, mesh->geometry(), active_cells, dofs, bs, fn,
560  coeffs, constant_values, cell_info);
561  }
562 
563  if (L.num_integrals(IntegralType::exterior_facet) > 0
564  or L.num_integrals(IntegralType::interior_facet) > 0)
565  {
566  // FIXME: cleanup these calls? Some of the happen internally again.
567  mesh->topology_mutable().create_entities(tdim - 1);
568  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
569  mesh->topology_mutable().create_entity_permutations();
570 
571  const std::vector<std::uint8_t>& perms
572  = mesh->topology().get_facet_permutations();
573  for (int i : L.integral_ids(IntegralType::exterior_facet))
574  {
575  const auto& fn = L.kernel(IntegralType::exterior_facet, i);
576  const std::vector<std::int32_t>& active_facets
577  = L.domains(IntegralType::exterior_facet, i);
578  fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, bs, fn,
579  coeffs, constant_values, cell_info,
580  perms);
581  }
582 
583  const std::vector<int> c_offsets = L.coefficient_offsets();
584  for (int i : L.integral_ids(IntegralType::interior_facet))
585  {
586  const auto& fn = L.kernel(IntegralType::interior_facet, i);
587  const std::vector<std::int32_t>& active_facets
588  = L.domains(IntegralType::interior_facet, i);
589  fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn,
590  coeffs, c_offsets, constant_values,
591  cell_info, perms);
592  }
593  }
594 }
595 //-----------------------------------------------------------------------------
596 template <typename T>
597 void assemble_cells(
598  tcb::span<T> b, const mesh::Geometry& geometry,
599  const std::vector<std::int32_t>& active_cells,
600  const graph::AdjacencyList<std::int32_t>& dofmap, const int bs,
601  const std::function<void(T*, const T*, const T*, const double*, const int*,
602  const std::uint8_t*, const std::uint32_t)>& kernel,
603  const array2d<T>& coeffs, const std::vector<T>& constant_values,
604  const std::vector<std::uint32_t>& cell_info)
605 {
606  const int gdim = geometry.dim();
607 
608  // Prepare cell geometry
609  const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
610 
611  // FIXME: Add proper interface for num coordinate dofs
612  const int num_dofs_g = x_dofmap.num_links(0);
613  const array2d<double>& x_g = geometry.x();
614 
615  // FIXME: Add proper interface for num_dofs
616  // Create data structures used in assembly
617  const int num_dofs = dofmap.links(0).size();
618  std::vector<double> coordinate_dofs(num_dofs_g * gdim);
619  std::vector<T> be(bs * num_dofs);
620 
621  // Iterate over active cells
622  for (std::int32_t c : active_cells)
623  {
624  // Get cell coordinates/geometry
625  auto x_dofs = x_dofmap.links(c);
626  for (std::size_t i = 0; i < x_dofs.size(); ++i)
627  {
628  std::copy_n(x_g.row(x_dofs[i]).begin(), gdim,
629  std::next(coordinate_dofs.begin(), i * gdim));
630  }
631 
632  // Tabulate vector for cell
633  std::fill(be.begin(), be.end(), 0);
634  kernel(be.data(), coeffs.row(c).data(), constant_values.data(),
635  coordinate_dofs.data(), nullptr, nullptr, cell_info[c]);
636 
637  // Scatter cell vector to 'global' vector array
638  auto dofs = dofmap.links(c);
639  for (int i = 0; i < num_dofs; ++i)
640  for (int k = 0; k < bs; ++k)
641  b[bs * dofs[i] + k] += be[bs * i + k];
642  }
643 }
644 //-----------------------------------------------------------------------------
645 template <typename T>
646 void assemble_exterior_facets(
647  tcb::span<T> b, const mesh::Mesh& mesh,
648  const std::vector<std::int32_t>& active_facets,
649  const graph::AdjacencyList<std::int32_t>& dofmap, const int bs,
650  const std::function<void(T*, const T*, const T*, const double*, const int*,
651  const std::uint8_t*, const std::uint32_t)>& fn,
652  const array2d<T>& coeffs, const std::vector<T>& constant_values,
653  const std::vector<std::uint32_t>& cell_info,
654  const std::vector<std::uint8_t>& perms)
655 {
656  const int gdim = mesh.geometry().dim();
657  const int tdim = mesh.topology().dim();
658 
659  // Prepare cell geometry
660  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
661 
662  // FIXME: Add proper interface for num coordinate dofs
663  const int num_dofs_g = x_dofmap.num_links(0);
664  const array2d<double>& x_g = mesh.geometry().x();
665 
666  // FIXME: Add proper interface for num_dofs
667  // Create data structures used in assembly
668  const int num_dofs = dofmap.links(0).size();
669  std::vector<double> coordinate_dofs(num_dofs_g * gdim);
670  std::vector<T> be(bs * num_dofs);
671 
672  auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
673  assert(f_to_c);
674  auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
675  assert(c_to_f);
676  for (std::int32_t f : active_facets)
677  {
678  // Get index of first attached cell
679  assert(f_to_c->num_links(f) > 0);
680  const std::int32_t cell = f_to_c->links(f)[0];
681 
682  // Get local index of facet with respect to the cell
683  auto facets = c_to_f->links(cell);
684  auto it = std::find(facets.begin(), facets.end(), f);
685  assert(it != facets.end());
686  const int local_facet = std::distance(facets.begin(), it);
687 
688  // Get cell coordinates/geometry
689  auto x_dofs = x_dofmap.links(cell);
690  for (std::size_t i = 0; i < x_dofs.size(); ++i)
691  {
692  std::copy_n(x_g.row(x_dofs[i]).begin(), gdim,
693  std::next(coordinate_dofs.begin(), i * gdim));
694  }
695 
696  // Tabulate element vector
697  std::fill(be.begin(), be.end(), 0);
698  fn(be.data(), coeffs.row(cell).data(), constant_values.data(),
699  coordinate_dofs.data(), &local_facet,
700  &perms[cell * facets.size() + local_facet], cell_info[cell]);
701 
702  // Add element vector to global vector
703  auto dofs = dofmap.links(cell);
704  for (int i = 0; i < num_dofs; ++i)
705  for (int k = 0; k < bs; ++k)
706  b[bs * dofs[i] + k] += be[bs * i + k];
707  }
708 }
709 //-----------------------------------------------------------------------------
710 template <typename T>
711 void assemble_interior_facets(
712  tcb::span<T> b, const mesh::Mesh& mesh,
713  const std::vector<std::int32_t>& active_facets, const fem::DofMap& dofmap,
714  const std::function<void(T*, const T*, const T*, const double*, const int*,
715  const std::uint8_t*, const std::uint32_t)>& fn,
716  const array2d<T>& coeffs, const std::vector<int>& offsets,
717  const std::vector<T>& constant_values,
718  const std::vector<std::uint32_t>& cell_info,
719  const std::vector<std::uint8_t>& perms)
720 {
721  const int gdim = mesh.geometry().dim();
722  const int tdim = mesh.topology().dim();
723 
724  // Prepare cell geometry
725  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
726 
727  // FIXME: Add proper interface for num coordinate dofs
728  const int num_dofs_g = x_dofmap.num_links(0);
729  const array2d<double>& x_g = mesh.geometry().x();
730 
731  // Create data structures used in assembly
732  std::vector<double> coordinate_dofs(2 * num_dofs_g * gdim);
733  std::vector<T> be;
734  std::vector<T> coeff_array(2 * offsets.back());
735  assert(offsets.back() == int(coeffs.shape[1]));
736 
737  const int bs = dofmap.bs();
738  auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
739  assert(f_to_c);
740  auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
741  assert(c_to_f);
742  const int offset_g = gdim * num_dofs_g;
743  for (std::int32_t f : active_facets)
744  {
745  // Get attached cell indices
746  auto cells = f_to_c->links(f);
747  assert(cells.size() == 2);
748 
749  const int facets_per_cell = c_to_f->num_links(cells[0]);
750 
751  // Create attached cells
752  std::array<int, 2> local_facet;
753  for (int i = 0; i < 2; ++i)
754  {
755  auto facets = c_to_f->links(cells[i]);
756  auto it = std::find(facets.begin(), facets.end(), f);
757  assert(it != facets.end());
758  local_facet[i] = std::distance(facets.begin(), it);
759  }
760 
761  // Get cell geometry
762  auto x_dofs0 = x_dofmap.links(cells[0]);
763  auto x_dofs1 = x_dofmap.links(cells[1]);
764  for (int i = 0; i < num_dofs_g; ++i)
765  {
766  for (int j = 0; j < gdim; ++j)
767  {
768  coordinate_dofs[i * gdim + j] = x_g(x_dofs0[i], j);
769  coordinate_dofs[offset_g + i * gdim + j] = x_g(x_dofs1[i], j);
770  }
771  }
772 
773  // Layout for the restricted coefficients is flattened
774  // w[coefficient][restriction][dof]
775  auto coeff_cell0 = coeffs.row(cells[0]);
776  auto coeff_cell1 = coeffs.row(cells[1]);
777 
778  // Loop over coefficients
779  for (std::size_t i = 0; i < offsets.size() - 1; ++i)
780  {
781  // Loop over entries for coefficient i
782  const int num_entries = offsets[i + 1] - offsets[i];
783  std::copy_n(coeff_cell0.data() + offsets[i], num_entries,
784  std::next(coeff_array.begin(), 2 * offsets[i]));
785  std::copy_n(coeff_cell1.data() + offsets[i], num_entries,
786  std::next(coeff_array.begin(), offsets[i + 1] + offsets[i]));
787  }
788 
789  // Get dofmaps for cells
790  tcb::span<const std::int32_t> dmap0 = dofmap.cell_dofs(cells[0]);
791  tcb::span<const std::int32_t> dmap1 = dofmap.cell_dofs(cells[1]);
792 
793  // Tabulate element vector
794  be.resize(bs * (dmap0.size() + dmap1.size()));
795  std::fill(be.begin(), be.end(), 0);
796  const std::array perm{perms[cells[0] * facets_per_cell + local_facet[0]],
797  perms[cells[1] * facets_per_cell + local_facet[1]]};
798  fn(be.data(), coeff_array.data(), constant_values.data(),
799  coordinate_dofs.data(), local_facet.data(), perm.data(),
800  cell_info[cells[0]]);
801 
802  // Add element vector to global vector
803  for (std::size_t i = 0; i < dmap0.size(); ++i)
804  for (int k = 0; k < bs; ++k)
805  b[bs * dmap0[i] + k] += be[bs * i + k];
806  for (std::size_t i = 0; i < dmap1.size(); ++i)
807  for (int k = 0; k < bs; ++k)
808  b[bs * dmap1[i] + k] += be[bs * (i + dmap0.size()) + k];
809  }
810 }
811 //-----------------------------------------------------------------------------
812 template <typename T>
813 void apply_lifting(
814  tcb::span<T> b, const std::vector<std::shared_ptr<const Form<T>>> a,
815  const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
816  const std::vector<tcb::span<const T>>& x0, double scale)
817 {
818  // FIXME: make changes to reactivate this check
819  if (!x0.empty() and x0.size() != a.size())
820  {
821  throw std::runtime_error(
822  "Mismatch in size between x0 and bilinear form in assembler.");
823  }
824 
825  if (a.size() != bcs1.size())
826  {
827  throw std::runtime_error(
828  "Mismatch in size between a and bcs in assembler.");
829  }
830 
831  for (std::size_t j = 0; j < a.size(); ++j)
832  {
833  std::vector<bool> bc_markers1;
834  std::vector<T> bc_values1;
835  if (a[j] and !bcs1[j].empty())
836  {
837  auto V1 = a[j]->function_spaces()[1];
838  assert(V1);
839  auto map1 = V1->dofmap()->index_map;
840  const int bs1 = V1->dofmap()->index_map_bs();
841  assert(map1);
842  const int crange = bs1 * (map1->size_local() + map1->num_ghosts());
843  bc_markers1.assign(crange, false);
844  bc_values1.assign(crange, 0.0);
845  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs1[j])
846  {
847  bc->mark_dofs(bc_markers1);
848  bc->dof_values(bc_values1);
849  }
850 
851  if (!x0.empty())
852  lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0[j], scale);
853  else
854  {
855  lift_bc<T>(b, *a[j], bc_values1, bc_markers1, tcb::span<const T>(),
856  scale);
857  }
858  }
859  }
860 }
861 //-----------------------------------------------------------------------------
862 template <typename T>
863 void lift_bc(tcb::span<T> b, const Form<T>& a,
864  const tcb::span<const T>& bc_values1,
865  const std::vector<bool>& bc_markers1, const tcb::span<const T>& x0,
866  double scale)
867 {
868  std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
869  assert(mesh);
870  const int tdim = mesh->topology().dim();
871  const std::int32_t num_cells
872  = mesh->topology().connectivity(tdim, 0)->num_nodes();
873 
874  // Get dofmap for columns and rows of a
875  assert(a.function_spaces().at(0));
876  assert(a.function_spaces().at(1));
877  const graph::AdjacencyList<std::int32_t>& dofmap0
878  = a.function_spaces()[0]->dofmap()->list();
879  const int bs0 = a.function_spaces()[0]->dofmap()->bs();
880  const graph::AdjacencyList<std::int32_t>& dofmap1
881  = a.function_spaces()[1]->dofmap()->list();
882  const int bs1 = a.function_spaces()[1]->dofmap()->bs();
883 
884  // Prepare constants
885  const std::vector<T> constant_values = pack_constants(a);
886 
887  // Prepare coefficients
888  const array2d<T> coeffs = pack_coefficients(a);
889 
890  const bool needs_permutation_data = a.needs_permutation_data();
891  if (needs_permutation_data)
892  mesh->topology_mutable().create_entity_permutations();
893  const std::vector<std::uint32_t>& cell_info
894  = needs_permutation_data ? mesh->topology().get_cell_permutation_info()
895  : std::vector<std::uint32_t>(num_cells);
896 
897  for (int i : a.integral_ids(IntegralType::cell))
898  {
899  const auto& kernel = a.kernel(IntegralType::cell, i);
900  const std::vector<std::int32_t>& active_cells
901  = a.domains(IntegralType::cell, i);
902  _lift_bc_cells(b, mesh->geometry(), kernel, active_cells, dofmap0, bs0,
903  dofmap1, bs1, coeffs, constant_values, cell_info, bc_values1,
904  bc_markers1, x0, scale);
905  }
906 
907  if (a.num_integrals(IntegralType::exterior_facet) > 0
908  or a.num_integrals(IntegralType::interior_facet) > 0)
909  {
910  // FIXME: cleanup these calls? Some of the happen internally again.
911  mesh->topology_mutable().create_entities(tdim - 1);
912  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
913  mesh->topology_mutable().create_entity_permutations();
914 
915  const std::vector<std::uint8_t>& perms
916  = mesh->topology().get_facet_permutations();
917  for (int i : a.integral_ids(IntegralType::exterior_facet))
918  {
919  const auto& kernel = a.kernel(IntegralType::exterior_facet, i);
920  const std::vector<std::int32_t>& active_facets
921  = a.domains(IntegralType::exterior_facet, i);
922  _lift_bc_exterior_facets(b, *mesh, kernel, active_facets, dofmap0, bs0,
923  dofmap1, bs1, coeffs, constant_values, cell_info,
924  perms, bc_values1, bc_markers1, x0, scale);
925  }
926 
927  const std::vector<int> c_offsets = a.coefficient_offsets();
928  for (int i : a.integral_ids(IntegralType::interior_facet))
929  {
930  const auto& kernel = a.kernel(IntegralType::interior_facet, i);
931  const std::vector<std::int32_t>& active_facets
932  = a.domains(IntegralType::interior_facet, i);
933  _lift_bc_interior_facets(b, *mesh, kernel, active_facets, dofmap0, bs0,
934  dofmap1, bs1, coeffs, c_offsets, constant_values,
935  cell_info, perms, bc_values1, bc_markers1, x0,
936  scale);
937  }
938  }
939 }
940 //-----------------------------------------------------------------------------
941 } // namespace dolfinx::fem::impl
void cells(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const std::reference_wrapper< const fem::DofMap >, 2 > &dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:18
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:453
void apply_lifting(tcb::span< T > 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< tcb::span< const T >> &x0, double scale)
Modify b such that:
Definition: assembler.h:69
array2d< typename U::scalar_type > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:398
void assemble_vector(tcb::span< T > b, const Form< T > &L)
Assemble linear form into a vector.
Definition: assembler.h:45