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