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