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