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 
117 template <typename T>
118 void lift_bc(
119  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const Form<T>& a,
120  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
121  const std::vector<bool>& bc_markers1, double scale);
122 
135 template <typename T>
136 void lift_bc(
137  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const Form<T>& a,
138  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
139  const std::vector<bool>& bc_markers1,
140  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
141  double scale);
142 
143 // Implementation of bc application
144 template <typename T>
145 void _lift_bc_cells(
146  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const Form<T>& a,
147  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
148  const std::vector<bool>& bc_markers1,
149  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
150  double scale)
151 {
152  assert(a.rank() == 2);
153 
154  // Get mesh from form
155  std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
156  assert(mesh);
157 
158  mesh->topology_mutable().create_entity_permutations();
159 
160  // Get dofmap for columns and rows of a
161  assert(a.function_space(0));
162  assert(a.function_space(1));
163  std::shared_ptr<const fem::DofMap> dofmap0 = a.function_space(0)->dofmap();
164  std::shared_ptr<const fem::DofMap> dofmap1 = a.function_space(1)->dofmap();
165  assert(dofmap0);
166  assert(dofmap1);
167 
168  // Prepare coefficients
169  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
170  = pack_coefficients(a);
171 
172  const std::function<void(T*, const T*, const T*, const double*, const int*,
173  const std::uint8_t*, const std::uint32_t)>& fn
174  = a.integrals().get_tabulate_tensor(IntegralType::cell, 0);
175 
176  // Prepare cell geometry
177  const int gdim = mesh->geometry().dim();
178  const graph::AdjacencyList<std::int32_t>& x_dofmap
179  = mesh->geometry().dofmap();
180 
181  // FIXME: Add proper interface for num coordinate dofs
182  const int num_dofs_g = x_dofmap.num_links(0);
183  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
184  = mesh->geometry().x();
185 
186  // Data structures used in bc application
187  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
188  coordinate_dofs(num_dofs_g, gdim);
189  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
190  Eigen::Matrix<T, Eigen::Dynamic, 1> be;
191 
192  // Prepare constants
193  if (!a.all_constants_set())
194  throw std::runtime_error("Unset constant in Form");
195  const Eigen::Array<T, Eigen::Dynamic, 1> constant_values = pack_constants(a);
196 
197  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
198  = mesh->topology().get_cell_permutation_info();
199 
200  // Iterate over all cells
201  const int tdim = mesh->topology().dim();
202  auto map = mesh->topology().index_map(tdim);
203  assert(map);
204  const int num_cells = map->size_local();
205  for (int c = 0; c < num_cells; ++c)
206  {
207  // Get dof maps for cell
208  auto dmap1 = dofmap1->cell_dofs(c);
209 
210  // Check if bc is applied to cell
211  bool has_bc = false;
212  for (Eigen::Index j = 0; j < dmap1.size(); ++j)
213  {
214  if (bc_markers1[dmap1[j]])
215  {
216  has_bc = true;
217  break;
218  }
219  }
220 
221  if (!has_bc)
222  continue;
223 
224  // Get cell vertex coordinates
225  auto x_dofs = x_dofmap.links(c);
226  for (int i = 0; i < num_dofs_g; ++i)
227  for (int j = 0; j < gdim; ++j)
228  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
229 
230  // Size data structure for assembly
231  auto dmap0 = dofmap0->cell_dofs(c);
232 
233  auto coeff_array = coeffs.row(c);
234  Ae.setZero(dmap0.size(), dmap1.size());
235  fn(Ae.data(), coeff_array.data(), constant_values.data(),
236  coordinate_dofs.data(), nullptr, nullptr, cell_info[c]);
237 
238  // Size data structure for assembly
239  be.setZero(dmap0.size());
240  for (Eigen::Index j = 0; j < dmap1.size(); ++j)
241  {
242  const std::int32_t jj = dmap1[j];
243  if (bc_markers1[jj])
244  {
245  const T bc = bc_values1[jj];
246  if (x0.rows() > 0)
247  be -= Ae.col(j) * scale * (bc - x0[jj]);
248  else
249  be -= Ae.col(j) * scale * bc;
250  }
251  }
252 
253  for (Eigen::Index k = 0; k < dmap0.size(); ++k)
254  b[dmap0[k]] += be[k];
255  }
256 }
257 //----------------------------------------------------------------------------
258 template <typename T>
259 void _lift_bc_exterior_facets(
260  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const Form<T>& a,
261  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
262  const std::vector<bool>& bc_markers1,
263  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
264  double scale)
265 {
266  assert(a.rank() == 2);
267 
268  // Get mesh from form
269  std::shared_ptr<const mesh::Mesh> mesh = a.mesh();
270  assert(mesh);
271 
272  mesh->topology_mutable().create_entity_permutations();
273 
274  const int gdim = mesh->geometry().dim();
275  const int tdim = mesh->topology().dim();
276 
277  // FIXME: cleanup these calls? Some of the happen internally again.
278  mesh->topology_mutable().create_entities(tdim - 1);
279  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
280  // FIXME: Why again -- appears already See 8 lines above.
281  mesh->topology_mutable().create_entity_permutations();
282 
283  // Get dofmap for columns and rows of a
284  assert(a.function_space(0));
285  assert(a.function_space(1));
286  std::shared_ptr<const fem::DofMap> dofmap0 = a.function_space(0)->dofmap();
287  std::shared_ptr<const fem::DofMap> dofmap1 = a.function_space(1)->dofmap();
288  assert(dofmap0);
289  assert(dofmap1);
290 
291  // Prepare coefficients
292  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
293  = pack_coefficients(a);
294 
295  const std::function<void(T*, const T*, const T*, const double*, const int*,
296  const std::uint8_t*, const std::uint32_t)>& fn
297  = a.integrals().get_tabulate_tensor(IntegralType::exterior_facet, 0);
298 
299  // Prepare cell geometry
300  const graph::AdjacencyList<std::int32_t>& x_dofmap
301  = mesh->geometry().dofmap();
302  // FIXME: Add proper interface for num coordinate dofs
303  const int num_dofs_g = x_dofmap.num_links(0);
304  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
305  = mesh->geometry().x();
306 
307  // Data structures used in bc application
308  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
309  coordinate_dofs(num_dofs_g, gdim);
310  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> Ae;
311  Eigen::Matrix<T, Eigen::Dynamic, 1> be;
312 
313  // Prepare constants
314  if (!a.all_constants_set())
315  throw std::runtime_error("Unset constant in Form");
316  const Eigen::Array<T, Eigen::Dynamic, 1> constant_values = pack_constants(a);
317 
318  // Iterate over owned facets
319  const mesh::Topology& topology = mesh->topology();
320  auto connectivity = topology.connectivity(tdim - 1, tdim);
321  assert(connectivity);
322  auto c_to_f = topology.connectivity(tdim, tdim - 1);
323  assert(c_to_f);
324  auto map = topology.index_map(tdim - 1);
325  assert(map);
326 
327  const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
328  = mesh->topology().get_facet_permutations();
329  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
330  = mesh->topology().get_cell_permutation_info();
331 
332  std::set<std::int32_t> fwd_shared_facets;
333  // Only need to consider shared facets when there are no ghost cells
334  if (topology.index_map(tdim)->num_ghosts() == 0)
335  {
336  fwd_shared_facets.insert(
337  topology.index_map(tdim - 1)->shared_indices().begin(),
338  topology.index_map(tdim - 1)->shared_indices().end());
339  }
340 
341  for (int f = 0; f < map->size_local(); ++f)
342  {
343  // Move to next facet if this one is an interior facet
344  // Interior facets have two attached cells. If on a process boundary,
345  // and owned locally, then they are "forward shared".
346  if (connectivity->num_links(f) == 2
347  or fwd_shared_facets.find(f) != fwd_shared_facets.end())
348  continue;
349 
350  // Create attached cell
351  assert(connectivity->num_links(f) == 1);
352  const std::int32_t cell = connectivity->links(f)[0];
353 
354  // Get local index of facet with respect to the cell
355  auto facets = c_to_f->links(cell);
356  const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
357  assert(it != (facets.data() + facets.rows()));
358  const int local_facet = std::distance(facets.data(), it);
359 
360  const std::uint8_t perm = perms(local_facet, cell);
361 
362  // Get dof maps for cell
363  auto dmap1 = dofmap1->cell_dofs(cell);
364 
365  // Check if bc is applied to cell
366  bool has_bc = false;
367  for (Eigen::Index j = 0; j < dmap1.size(); ++j)
368  {
369  if (bc_markers1[dmap1[j]])
370  {
371  has_bc = true;
372  break;
373  }
374  }
375 
376  if (!has_bc)
377  continue;
378 
379  // Get cell vertex coordinates
380  auto x_dofs = x_dofmap.links(cell);
381  for (int i = 0; i < num_dofs_g; ++i)
382  for (int j = 0; j < gdim; ++j)
383  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
384 
385  // Size data structure for assembly
386  auto dmap0 = dofmap0->cell_dofs(cell);
387 
388  // TODO: Move gathering of coefficients outside of main assembly
389  // loop
390 
391  auto coeff_array = coeffs.row(cell);
392  Ae.setZero(dmap0.size(), dmap1.size());
393  fn(Ae.data(), coeff_array.data(), constant_values.data(),
394  coordinate_dofs.data(), &local_facet, &perm, cell_info[cell]);
395 
396  // Size data structure for assembly
397  be.setZero(dmap0.size());
398  for (Eigen::Index j = 0; j < dmap1.size(); ++j)
399  {
400  const std::int32_t jj = dmap1[j];
401  if (bc_markers1[jj])
402  {
403  const T bc = bc_values1[jj];
404  if (x0.rows() > 0)
405  be -= Ae.col(j) * scale * (bc - x0[jj]);
406  else
407  be -= Ae.col(j) * scale * bc;
408  }
409  }
410 
411  for (Eigen::Index k = 0; k < dmap0.size(); ++k)
412  b[dmap0[k]] += be[k];
413  }
414 }
415 
416 //-----------------------------------------------------------------------------
417 template <typename T>
418 void assemble_vector(Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
419  const Form<T>& L)
420 {
421  std::shared_ptr<const mesh::Mesh> mesh = L.mesh();
422  assert(mesh);
423  const int tdim = mesh->topology().dim();
424  const std::int32_t num_cells
425  = mesh->topology().connectivity(tdim, 0)->num_nodes();
426 
427  // Get dofmap data
428  assert(L.function_space(0));
429  std::shared_ptr<const fem::DofMap> dofmap = L.function_space(0)->dofmap();
430  assert(dofmap);
431  const graph::AdjacencyList<std::int32_t>& dofs = dofmap->list();
432 
433  // Prepare constants
434  if (!L.all_constants_set())
435  throw std::runtime_error("Unset constant in Form");
436  const Eigen::Array<T, Eigen::Dynamic, 1> constant_values = pack_constants(L);
437 
438  // Prepare coefficients
439  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
440  = pack_coefficients(L);
441 
442  const FormIntegrals<T>& integrals = L.integrals();
443  const bool needs_permutation_data = integrals.needs_permutation_data();
444  if (needs_permutation_data)
445  mesh->topology_mutable().create_entity_permutations();
446  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
447  = needs_permutation_data
448  ? mesh->topology().get_cell_permutation_info()
449  : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
450 
451  for (int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i)
452  {
453  const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i);
454  const std::vector<std::int32_t>& active_cells
455  = integrals.integral_domains(IntegralType::cell, i);
456  fem::impl::assemble_cells(b, mesh->geometry(), active_cells, dofs, fn,
457  coeffs, constant_values, cell_info);
458  }
459 
460  if (integrals.num_integrals(IntegralType::exterior_facet) > 0
461  or integrals.num_integrals(IntegralType::interior_facet) > 0)
462  {
463  // FIXME: cleanup these calls? Some of the happen internally again.
464  mesh->topology_mutable().create_entities(tdim - 1);
465  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
466 
467  const int facets_per_cell
468  = mesh::cell_num_entities(mesh->topology().cell_type(), tdim - 1);
469  const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
470  = needs_permutation_data
471  ? mesh->topology().get_facet_permutations()
472  : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
473  facets_per_cell, num_cells);
474  for (int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet);
475  ++i)
476  {
477  const auto& fn
478  = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i);
479  const std::vector<std::int32_t>& active_facets
480  = integrals.integral_domains(IntegralType::exterior_facet, i);
481  fem::impl::assemble_exterior_facets(b, *mesh, active_facets, dofs, fn,
482  coeffs, constant_values, cell_info,
483  perms);
484  }
485 
486  const std::vector<int> c_offsets = L.coefficients().offsets();
487  for (int i = 0; i < integrals.num_integrals(IntegralType::interior_facet);
488  ++i)
489  {
490  const auto& fn
491  = integrals.get_tabulate_tensor(IntegralType::interior_facet, i);
492  const std::vector<std::int32_t>& active_facets
493  = integrals.integral_domains(IntegralType::interior_facet, i);
494  fem::impl::assemble_interior_facets(b, *mesh, active_facets, *dofmap, fn,
495  coeffs, c_offsets, constant_values,
496  cell_info, perms);
497  }
498  }
499 }
500 //-----------------------------------------------------------------------------
501 template <typename T>
502 void assemble_cells(
503  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
504  const mesh::Geometry& geometry,
505  const std::vector<std::int32_t>& active_cells,
506  const graph::AdjacencyList<std::int32_t>& dofmap,
507  const std::function<void(T*, const T*, const T*, const double*, const int*,
508  const std::uint8_t*, const std::uint32_t)>& kernel,
509  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
510  coeffs,
511  const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
512  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info)
513 {
514  const int gdim = geometry.dim();
515 
516  // Prepare cell geometry
517  const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
518 
519  // FIXME: Add proper interface for num coordinate dofs
520  const int num_dofs_g = x_dofmap.num_links(0);
521  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
522  = geometry.x();
523 
524  // FIXME: Add proper interface for num_dofs
525  // Create data structures used in assembly
526  const int num_dofs = dofmap.links(0).size();
527  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
528  coordinate_dofs(num_dofs_g, gdim);
529  Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
530 
531  // Iterate over active cells
532  for (std::int32_t c : active_cells)
533  {
534  // Get cell coordinates/geometry
535  auto x_dofs = x_dofmap.links(c);
536  for (int i = 0; i < num_dofs_g; ++i)
537  for (int j = 0; j < gdim; ++j)
538  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
539 
540  // Tabulate vector for cell
541  std::fill(be.data(), be.data() + num_dofs, 0);
542  kernel(be.data(), coeffs.row(c).data(), constant_values.data(),
543  coordinate_dofs.data(), nullptr, nullptr, cell_info[c]);
544 
545  // Scatter cell vector to 'global' vector array
546  auto dofs = dofmap.links(c);
547  for (Eigen::Index i = 0; i < num_dofs; ++i)
548  b[dofs[i]] += be[i];
549  }
550 }
551 //-----------------------------------------------------------------------------
552 template <typename T>
553 void assemble_exterior_facets(
554  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const mesh::Mesh& mesh,
555  const std::vector<std::int32_t>& active_facets,
556  const graph::AdjacencyList<std::int32_t>& dofmap,
557  const std::function<void(T*, const T*, const T*, const double*, const int*,
558  const std::uint8_t*, const std::uint32_t)>& fn,
559  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
560  coeffs,
561  const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
562  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
563  const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
564 {
565  const int gdim = mesh.geometry().dim();
566  const int tdim = mesh.topology().dim();
567 
568  // Prepare cell geometry
569  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
570 
571  // FIXME: Add proper interface for num coordinate dofs
572  const int num_dofs_g = x_dofmap.num_links(0);
573  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
574  = mesh.geometry().x();
575 
576  // FIXME: Add proper interface for num_dofs
577  // Create data structures used in assembly
578  const int num_dofs = dofmap.links(0).size();
579  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
580  coordinate_dofs(num_dofs_g, gdim);
581  Eigen::Matrix<T, Eigen::Dynamic, 1> be(num_dofs);
582 
583  auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
584  assert(f_to_c);
585  auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
586  assert(c_to_f);
587  for (const auto& f : active_facets)
588  {
589  // Get index of first attached cell
590  assert(f_to_c->num_links(f) > 0);
591  const std::int32_t cell = f_to_c->links(f)[0];
592 
593  // Get local index of facet with respect to the cell
594  auto facets = c_to_f->links(cell);
595  const auto* it = std::find(facets.data(), facets.data() + facets.rows(), f);
596  assert(it != (facets.data() + facets.rows()));
597  const int local_facet = std::distance(facets.data(), it);
598 
599  // Get cell coordinates/geometry
600  auto x_dofs = x_dofmap.links(cell);
601  for (int i = 0; i < num_dofs_g; ++i)
602  for (int j = 0; j < gdim; ++j)
603  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
604 
605  // Tabulate element vector
606  std::fill(be.data(), be.data() + num_dofs, 0);
607  fn(be.data(), coeffs.row(cell).data(), constant_values.data(),
608  coordinate_dofs.data(), &local_facet, &perms(local_facet, cell),
609  cell_info[cell]);
610 
611  // Add element vector to global vector
612  auto dofs = dofmap.links(cell);
613  for (Eigen::Index i = 0; i < num_dofs; ++i)
614  b[dofs[i]] += be[i];
615  }
616 }
617 //-----------------------------------------------------------------------------
618 template <typename T>
619 void assemble_interior_facets(
620  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const mesh::Mesh& mesh,
621  const std::vector<std::int32_t>& active_facets, const fem::DofMap& dofmap,
622  const std::function<void(T*, const T*, const T*, const double*, const int*,
623  const std::uint8_t*, const std::uint32_t)>& fn,
624  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
625  coeffs,
626  const std::vector<int>& offsets,
627  const Eigen::Array<T, Eigen::Dynamic, 1>& constant_values,
628  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
629  const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
630 {
631  const int gdim = mesh.geometry().dim();
632  const int tdim = mesh.topology().dim();
633 
634  // Prepare cell geometry
635  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
636  // FIXME: Add proper interface for num coordinate dofs
637 
638  const int num_dofs_g = x_dofmap.num_links(0);
639  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
640  = mesh.geometry().x();
641 
642  // Creat data structures used in assembly
643  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
644  coordinate_dofs(2 * num_dofs_g, gdim);
645  Eigen::Matrix<T, Eigen::Dynamic, 1> be;
646  Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
647  assert(offsets.back() == coeffs.cols());
648 
649  auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
650  assert(f_to_c);
651  auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
652  assert(c_to_f);
653  for (const auto& f : active_facets)
654  {
655  // Get attached cell indices
656  auto cells = f_to_c->links(f);
657  assert(cells.rows() == 2);
658 
659  // Create attached cells
660  std::array<int, 2> local_facet;
661  for (int i = 0; i < 2; ++i)
662  {
663  auto facets = c_to_f->links(cells[i]);
664  const auto* it
665  = std::find(facets.data(), facets.data() + facets.rows(), f);
666  assert(it != (facets.data() + facets.rows()));
667  local_facet[i] = std::distance(facets.data(), it);
668  }
669 
670  // Get cell geometry
671  auto x_dofs0 = x_dofmap.links(cells[0]);
672  auto x_dofs1 = x_dofmap.links(cells[1]);
673  for (int i = 0; i < num_dofs_g; ++i)
674  {
675  for (int j = 0; j < gdim; ++j)
676  {
677  coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
678  coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
679  }
680  }
681 
682  // Get dofmaps for cell
683  auto dmap0 = dofmap.cell_dofs(cells[0]);
684  auto dmap1 = dofmap.cell_dofs(cells[1]);
685 
686  // Layout for the restricted coefficients is flattened
687  // w[coefficient][restriction][dof]
688  auto coeff_cell0 = coeffs.row(cells[0]);
689  auto coeff_cell1 = coeffs.row(cells[1]);
690 
691  // Loop over coefficients
692  for (std::size_t i = 0; i < offsets.size() - 1; ++i)
693  {
694  // Loop over entries for coefficient i
695  const int num_entries = offsets[i + 1] - offsets[i];
696  coeff_array.segment(2 * offsets[i], num_entries)
697  = coeff_cell0.segment(offsets[i], num_entries);
698  coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
699  = coeff_cell1.segment(offsets[i], num_entries);
700  }
701 
702  // Tabulate element vector
703  be.setZero(dmap0.size() + dmap1.size());
704 
705  const std::array perm{perms(local_facet[0], cells[0]),
706  perms(local_facet[1], cells[1])};
707  fn(be.data(), coeff_array.data(), constant_values.data(),
708  coordinate_dofs.data(), local_facet.data(), perm.data(),
709  cell_info[cells[0]]);
710 
711  // Add element vector to global vector
712  for (Eigen::Index i = 0; i < dmap0.size(); ++i)
713  b[dmap0[i]] += be[i];
714  for (Eigen::Index i = 0; i < dmap1.size(); ++i)
715  b[dmap1[i]] += be[i + dmap0.size()];
716  }
717 }
718 //-----------------------------------------------------------------------------
719 template <typename T>
720 void apply_lifting(
721  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b,
722  const std::vector<std::shared_ptr<const Form<T>>> a,
723  const std::vector<std::vector<std::shared_ptr<const DirichletBC<T>>>>& bcs1,
724  const std::vector<Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>>&
725  x0,
726  double scale)
727 {
728  // FIXME: make changes to reactivate this check
729  if (!x0.empty() and x0.size() != a.size())
730  {
731  throw std::runtime_error(
732  "Mismatch in size between x0 and bilinear form in assembler.");
733  }
734 
735  if (a.size() != bcs1.size())
736  {
737  throw std::runtime_error(
738  "Mismatch in size between a and bcs in assembler.");
739  }
740 
741  for (std::size_t j = 0; j < a.size(); ++j)
742  {
743  std::vector<bool> bc_markers1;
744  Eigen::Matrix<T, Eigen::Dynamic, 1> bc_values1;
745  if (a[j] and !bcs1[j].empty())
746  {
747  auto V1 = a[j]->function_space(1);
748  assert(V1);
749  auto map1 = V1->dofmap()->index_map;
750  assert(map1);
751  const int crange
752  = map1->block_size() * (map1->size_local() + map1->num_ghosts());
753  bc_markers1.assign(crange, false);
754  bc_values1 = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(crange);
755  for (const std::shared_ptr<const DirichletBC<T>>& bc : bcs1[j])
756  {
757  bc->mark_dofs(bc_markers1);
758  bc->dof_values(bc_values1);
759  }
760 
761  // Modify (apply lifting) vector
762  if (!x0.empty())
763  lift_bc<T>(b, *a[j], bc_values1, bc_markers1, x0[j], scale);
764  else
765  lift_bc<T>(b, *a[j], bc_values1, bc_markers1, scale);
766  }
767  }
768 }
769 //-----------------------------------------------------------------------------
770 template <typename T>
771 void lift_bc(
772  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const Form<T>& a,
773  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
774  const std::vector<bool>& bc_markers1, double scale)
775 {
776  const Eigen::Matrix<T, Eigen::Dynamic, 1> x0(0);
777  if (a.integrals().num_integrals(fem::IntegralType::cell) > 0)
778  _lift_bc_cells<T>(b, a, bc_values1, bc_markers1, x0, scale);
779  if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0)
780  _lift_bc_exterior_facets<T>(b, a, bc_values1, bc_markers1, x0, scale);
781 }
782 //-----------------------------------------------------------------------------
783 template <typename T>
784 void lift_bc(
785  Eigen::Ref<Eigen::Matrix<T, Eigen::Dynamic, 1>> b, const Form<T>& a,
786  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& bc_values1,
787  const std::vector<bool>& bc_markers1,
788  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& x0,
789  double scale)
790 {
791  if (a.integrals().num_integrals(fem::IntegralType::cell) > 0)
792  _lift_bc_cells(b, a, bc_values1, bc_markers1, x0, scale);
793  if (a.integrals().num_integrals(fem::IntegralType::exterior_facet) > 0)
794  _lift_bc_exterior_facets(b, a, bc_values1, bc_markers1, x0, scale);
795 }
796 //-----------------------------------------------------------------------------
797 } // namespace dolfinx::fem::impl
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::pack_constants
Eigen::Array< T, Eigen::Dynamic, 1 > pack_constants(const fem::Form< T > &form)
Pack form constants ready for assembly.
Definition: utils.h:359
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_coefficients
Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > pack_coefficients(const fem::Form< T > &form)
Pack form coefficients ready for assembly.
Definition: utils.h:314