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