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