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