DOLFIN-X
DOLFIN-X C++ interface
assemble_scalar_impl.h
1 // Copyright (C) 2019-2020 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 "Form.h"
10 #include "utils.h"
11 #include <Eigen/Dense>
12 #include <dolfinx/common/IndexMap.h>
13 #include <dolfinx/common/types.h>
14 #include <dolfinx/function/Constant.h>
15 #include <dolfinx/function/FunctionSpace.h>
16 #include <dolfinx/mesh/Geometry.h>
17 #include <dolfinx/mesh/Mesh.h>
18 #include <dolfinx/mesh/Topology.h>
19 #include <memory>
20 #include <vector>
21 
22 namespace dolfinx::fem::impl
23 {
24 
26 template <typename T>
27 T assemble_scalar(const fem::Form<T>& M);
28 
30 template <typename T>
31 T assemble_cells(
32  const mesh::Mesh& mesh, const std::vector<std::int32_t>& active_cells,
33  const std::function<void(T*, const T*, const T*, const double*, const int*,
34  const std::uint8_t*, const std::uint32_t)>& fn,
35  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
36  coeffs,
37  const std::vector<T>& constant_values);
38 
40 template <typename T>
41 T assemble_exterior_facets(
42  const mesh::Mesh& mesh, const std::vector<std::int32_t>& active_cells,
43  const std::function<void(T*, const T*, const T*, const double*, const int*,
44  const std::uint8_t*, const std::uint32_t)>& fn,
45  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
46  coeffs,
47  const std::vector<T>& constant_values);
48 
50 template <typename T>
51 T assemble_interior_facets(
52  const mesh::Mesh& mesh, const std::vector<std::int32_t>& active_cells,
53  const std::function<void(T*, const T*, const T*, const double*, const int*,
54  const std::uint8_t*, const std::uint32_t)>& fn,
55  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
56  coeffs,
57  const std::vector<int>& offsets, const std::vector<T>& constant_values);
58 
59 //-----------------------------------------------------------------------------
60 template <typename T>
61 T assemble_scalar(const fem::Form<T>& M)
62 {
63  std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
64  assert(mesh);
65 
66  // Prepare constants
67  if (!M.all_constants_set())
68  throw std::runtime_error("Unset constant in Form");
69  auto constants = M.constants();
70 
71  std::vector<T> constant_values;
72  for (auto const& constant : constants)
73  {
74  // Get underlying data array of this Constant
75  const std::vector<T>& array = constant.second->value;
76  constant_values.insert(constant_values.end(), array.data(),
77  array.data() + array.size());
78  }
79 
80  // Prepare coefficients
81  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
82  = pack_coefficients(M);
83 
84  const FormIntegrals<T>& integrals = M.integrals();
85  T value(0);
86  for (int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i)
87  {
88  const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i);
89  const std::vector<std::int32_t>& active_cells
90  = integrals.integral_domains(IntegralType::cell, i);
91  value += fem::impl::assemble_cells(*mesh, active_cells, fn, coeffs,
92  constant_values);
93  }
94 
95  for (int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet);
96  ++i)
97  {
98  const auto& fn
99  = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i);
100  const std::vector<std::int32_t>& active_facets
101  = integrals.integral_domains(IntegralType::exterior_facet, i);
102  value += fem::impl::assemble_exterior_facets(*mesh, active_facets, fn,
103  coeffs, constant_values);
104  }
105 
106  const std::vector<int> c_offsets = M.coefficients().offsets();
107  for (int i = 0; i < integrals.num_integrals(IntegralType::interior_facet);
108  ++i)
109  {
110  const auto& fn
111  = integrals.get_tabulate_tensor(IntegralType::interior_facet, i);
112  const std::vector<std::int32_t>& active_facets
113  = integrals.integral_domains(IntegralType::interior_facet, i);
114  value += fem::impl::assemble_interior_facets(
115  *mesh, active_facets, fn, coeffs, c_offsets, constant_values);
116  }
117 
118  return value;
119 }
120 //-----------------------------------------------------------------------------
121 template <typename T>
122 T assemble_cells(
123  const mesh::Mesh& mesh, const std::vector<std::int32_t>& active_cells,
124  const std::function<void(T*, const T*, const T*, const double*, const int*,
125  const std::uint8_t*, const std::uint32_t)>& fn,
126  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
127  coeffs,
128  const std::vector<T>& constant_values)
129 {
130  const int gdim = mesh.geometry().dim();
131  const int tdim = mesh.topology().dim();
132  mesh.topology_mutable().create_entities(tdim);
133  mesh.topology_mutable().create_entity_permutations();
134 
135  // Prepare cell geometry
136  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
137 
138  // FIXME: Add proper interface for num coordinate dofs
139  const int num_dofs_g = x_dofmap.num_links(0);
140  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
141  = mesh.geometry().x();
142 
143  // Create data structures used in assembly
144  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
145  coordinate_dofs(num_dofs_g, gdim);
146 
147  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
148  = mesh.topology().get_cell_permutation_info();
149 
150  // Iterate over all cells
151  T value(0);
152  for (std::int32_t c : active_cells)
153  {
154  // Get cell coordinates/geometry
155  auto x_dofs = x_dofmap.links(c);
156  for (int i = 0; i < num_dofs_g; ++i)
157  for (int j = 0; j < gdim; ++j)
158  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
159 
160  auto coeff_cell = coeffs.row(c);
161  fn(&value, coeff_cell.data(), constant_values.data(),
162  coordinate_dofs.data(), nullptr, nullptr, cell_info[c]);
163  }
164 
165  return value;
166 }
167 //-----------------------------------------------------------------------------
168 template <typename T>
169 T assemble_exterior_facets(
170  const mesh::Mesh& mesh, const std::vector<std::int32_t>& active_facets,
171  const std::function<void(T*, const T*, const T*, const double*, const int*,
172  const std::uint8_t*, const std::uint32_t)>& fn,
173  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
174  coeffs,
175  const std::vector<T>& constant_values)
176 {
177  const int gdim = mesh.geometry().dim();
178  const int tdim = mesh.topology().dim();
179 
180  // FIXME: cleanup these calls? Some of these happen internally again.
181  mesh.topology_mutable().create_entities(tdim - 1);
182  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
183  mesh.topology_mutable().create_entity_permutations();
184 
185  // Prepare cell geometry
186  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
187 
188  // FIXME: Add proper interface for num coordinate dofs
189  const int num_dofs_g = x_dofmap.num_links(0);
190  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
191  = mesh.geometry().x();
192 
193  // Creat data structures used in assembly
194  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
195  coordinate_dofs(num_dofs_g, gdim);
196 
197  const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
198  = mesh.topology().get_facet_permutations();
199  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
200  = mesh.topology().get_cell_permutation_info();
201 
202  auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
203  assert(f_to_c);
204  auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
205  assert(c_to_f);
206 
207  // Iterate over all facets
208  T value(0);
209  for (std::int32_t facet : active_facets)
210  {
211  // Create attached cell
212  assert(f_to_c->num_links(facet) == 1);
213  const int cell = f_to_c->links(facet)[0];
214 
215  // Get local index of facet with respect to the cell
216  auto facets = c_to_f->links(cell);
217  const auto* it
218  = std::find(facets.data(), facets.data() + facets.rows(), facet);
219  assert(it != (facets.data() + facets.rows()));
220  const int local_facet = std::distance(facets.data(), it);
221 
222  auto x_dofs = x_dofmap.links(cell);
223  for (int i = 0; i < num_dofs_g; ++i)
224  for (int j = 0; j < gdim; ++j)
225  coordinate_dofs(i, j) = x_g(x_dofs[i], j);
226 
227  auto coeff_cell = coeffs.row(cell);
228  const std::uint8_t perm = perms(local_facet, cell);
229  fn(&value, coeff_cell.data(), constant_values.data(),
230  coordinate_dofs.data(), &local_facet, &perm, cell_info[cell]);
231  }
232 
233  return value;
234 }
235 //-----------------------------------------------------------------------------
236 template <typename T>
237 T assemble_interior_facets(
238  const mesh::Mesh& mesh, const std::vector<std::int32_t>& active_facets,
239  const std::function<void(T*, const T*, const T*, const double*, const int*,
240  const std::uint8_t*, const std::uint32_t)>& fn,
241  const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
242  coeffs,
243  const std::vector<int>& offsets, const std::vector<T>& constant_values)
244 {
245  const int gdim = mesh.geometry().dim();
246  const int tdim = mesh.topology().dim();
247 
248  // FIXME: cleanup these calls? Some of the happen internally again.
249  mesh.topology_mutable().create_entities(tdim - 1);
250  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
251  mesh.topology_mutable().create_entity_permutations();
252 
253  // Prepare cell geometry
254  const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
255 
256  // FIXME: Add proper interface for num coordinate dofs
257  const int num_dofs_g = x_dofmap.num_links(0);
258  const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
259  = mesh.geometry().x();
260 
261  // Creat data structures used in assembly
262  Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
263  coordinate_dofs(2 * num_dofs_g, gdim);
264  Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
265  assert(offsets.back() == coeffs.cols());
266 
267  const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
268  = mesh.topology().get_facet_permutations();
269  const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
270  = mesh.topology().get_cell_permutation_info();
271 
272  auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
273  assert(f_to_c);
274  auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
275  assert(c_to_f);
276 
277  // Iterate over all facets
278  T value(0);
279  for (std::int32_t f : active_facets)
280  {
281  // Create attached cell
282  auto cells = f_to_c->links(f);
283  assert(cells.rows() == 2);
284 
285  // Get local index of facet with respect to the cell
286  std::array<int, 2> local_facet;
287  for (int i = 0; i < 2; ++i)
288  {
289  auto facets = c_to_f->links(cells[i]);
290  const auto* it
291  = std::find(facets.data(), facets.data() + facets.rows(), f);
292  assert(it != (facets.data() + facets.rows()));
293  local_facet[i] = std::distance(facets.data(), it);
294  }
295 
296  const std::array perm{perms(local_facet[0], cells[0]),
297  perms(local_facet[1], cells[1])};
298 
299  // Get cell geometry
300  auto x_dofs0 = x_dofmap.links(cells[0]);
301  auto x_dofs1 = x_dofmap.links(cells[1]);
302  for (int i = 0; i < num_dofs_g; ++i)
303  {
304  for (int j = 0; j < gdim; ++j)
305  {
306  coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
307  coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
308  }
309  }
310 
311  // Layout for the restricted coefficients is flattened
312  // w[coefficient][restriction][dof]
313  auto coeff_cell0 = coeffs.row(cells[0]);
314  auto coeff_cell1 = coeffs.row(cells[1]);
315 
316  // Loop over coefficients
317  for (std::size_t i = 0; i < offsets.size() - 1; ++i)
318  {
319  // Loop over entries for coefficient i
320  const int num_entries = offsets[i + 1] - offsets[i];
321  coeff_array.segment(2 * offsets[i], num_entries)
322  = coeff_cell0.segment(offsets[i], num_entries);
323  coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
324  = coeff_cell1.segment(offsets[i], num_entries);
325  }
326  fn(&value, coeff_array.data(), constant_values.data(),
327  coordinate_dofs.data(), local_facet.data(), perm.data(),
328  cell_info[cells[0]]);
329  }
330 
331  return value;
332 }
333 //-----------------------------------------------------------------------------
334 
335 } // namespace dolfinx::fem::impl
dolfinx::fem::assemble_scalar
T assemble_scalar(const Form< T > &M)
Assemble functional into scalar. Caller is responsible for accumulation across processes.
Definition: assembler.h:39
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