DOLFIN-X
DOLFIN-X C++ interface
FormIntegrals.h
1 // Copyright (C) 2018 Chris Richardson
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 <array>
10 #include <dolfinx/mesh/MeshTags.h>
11 #include <functional>
12 #include <set>
13 #include <vector>
14 
15 namespace dolfinx
16 {
17 namespace mesh
18 {
19 class Mesh;
20 } // namespace mesh
21 
22 namespace fem
23 {
24 
26 enum class IntegralType : std::int8_t
27 {
28  cell = 0,
29  exterior_facet = 1,
30  interior_facet = 2,
31  vertex = 3
32 };
33 
36 
37 template <typename T>
39 {
40 public:
49  const std::map<
51  std::vector<std::pair<
52  int, std::function<void(T*, const T*, const T*, const double*,
53  const int*, const std::uint8_t*,
54  const std::uint32_t)>>>>& integrals,
56  : _needs_permutation_data(needs_permutation_data)
57  {
58  for (auto& integral_type : integrals)
59  {
60  for (auto& integral : integral_type.second)
61  {
62  set_tabulate_tensor(integral_type.first, integral.first,
63  integral.second);
64  }
65  }
66  };
67 
73  const std::function<void(T*, const T*, const T*, const double*, const int*,
74  const std::uint8_t*, const std::uint32_t)>&
75  get_tabulate_tensor(IntegralType type, int i) const
76  {
77  return _integrals.at(static_cast<int>(type)).at(i).tabulate;
78  }
79 
88  IntegralType type, int i,
89  const std::function<void(T*, const T*, const T*, const double*,
90  const int*, const std::uint8_t*,
91  const std::uint32_t)>& fn)
92  {
93  std::vector<struct FormIntegrals::Integral>& integrals
94  = _integrals.at(static_cast<int>(type));
95 
96  // Find insertion point
97  int pos = 0;
98  for (const auto& q : integrals)
99  {
100  if (q.id == i)
101  {
102  throw std::runtime_error("Integral with ID " + std::to_string(i)
103  + " already exists");
104  }
105  else if (q.id > i)
106  break;
107  ++pos;
108  }
109 
110  // Insert new Integral
111  integrals.insert(integrals.begin() + pos, {fn, i, {}});
112  }
113 
116  std::set<IntegralType> types() const
117  {
118  static const std::array types{IntegralType::cell,
119  IntegralType::exterior_facet,
120  IntegralType::interior_facet};
121  std::set<IntegralType> set;
122  for (auto type : types)
123  if (!_integrals.at(static_cast<int>(type)).empty())
124  set.insert(type);
125  return set;
126  }
127 
131  int num_integrals(IntegralType type) const
132  {
133  return _integrals.at(static_cast<int>(type)).size();
134  }
135 
141  std::vector<int> integral_ids(IntegralType type) const
142  {
143  std::vector<int> ids;
144  for (const auto& integral : _integrals[static_cast<int>(type)])
145  ids.push_back(integral.id);
146  return ids;
147  }
148 
157  const std::vector<std::int32_t>& integral_domains(IntegralType type,
158  int i) const
159  {
160  return _integrals.at(static_cast<int>(type)).at(i).active_entities;
161  }
162 
168  void set_domains(IntegralType type, const mesh::MeshTags<int>& marker)
169  {
170  std::vector<struct Integral>& integrals
171  = _integrals.at(static_cast<int>(type));
172  if (integrals.size() == 0)
173  return;
174 
175  std::shared_ptr<const mesh::Mesh> mesh = marker.mesh();
176  const mesh::Topology& topology = mesh->topology();
177  const int tdim = topology.dim();
178  int dim = tdim;
179  if (type == IntegralType::exterior_facet
180  or type == IntegralType::interior_facet)
181  {
182  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
183  dim = tdim - 1;
184  }
185  else if (type == IntegralType::vertex)
186  dim = 0;
187 
188  if (dim != marker.dim())
189  {
190  throw std::runtime_error("Invalid MeshTags dimension:"
191  + std::to_string(marker.dim()));
192  }
193 
194  // Create a reverse map
195  std::map<int, int> id_to_integral;
196  for (std::size_t i = 0; i < integrals.size(); ++i)
197  {
198  if (integrals[i].id != -1)
199  {
200  integrals[i].active_entities.clear();
201  id_to_integral.insert({integrals[i].id, i});
202  }
203  }
204 
205  // Get mesh tag data
206  const std::vector<int>& values = marker.values();
207  const std::vector<std::int32_t>& tagged_entities = marker.indices();
208  assert(topology.index_map(dim));
209  const auto entity_end
210  = std::lower_bound(tagged_entities.begin(), tagged_entities.end(),
211  topology.index_map(dim)->size_local());
212 
213  if (type == IntegralType::exterior_facet)
214  {
215  auto f_to_c = topology.connectivity(tdim - 1, tdim);
216  assert(f_to_c);
217  if (type == IntegralType::exterior_facet)
218  {
219  // Only need to consider shared facets when there are no ghost cells
220  assert(topology.index_map(tdim));
221  std::set<std::int32_t> fwd_shared;
222  if (topology.index_map(tdim)->num_ghosts() == 0)
223  {
224  fwd_shared.insert(
225  topology.index_map(tdim - 1)->shared_indices().begin(),
226  topology.index_map(tdim - 1)->shared_indices().end());
227  }
228 
229  for (auto f = tagged_entities.begin(); f != entity_end; ++f)
230  {
231  const std::size_t i = std::distance(tagged_entities.cbegin(), f);
232 
233  // All "owned" facets connected to one cell, that are not shared,
234  // should be external.
235  if (f_to_c->num_links(*f) == 1
236  and fwd_shared.find(*f) == fwd_shared.end())
237  {
238  if (auto it = id_to_integral.find(values[i]);
239  it != id_to_integral.end())
240  {
241  integrals[it->second].active_entities.push_back(*f);
242  }
243  }
244  }
245  }
246  else if (type == IntegralType::interior_facet)
247  {
248  for (auto f = tagged_entities.begin(); f != entity_end; ++f)
249  {
250  if (f_to_c->num_links(*f) == 2)
251  {
252  const std::size_t i = std::distance(tagged_entities.cbegin(), f);
253  if (const auto it = id_to_integral.find(values[i]);
254  it != id_to_integral.end())
255  {
256  integrals[it->second].active_entities.push_back(*f);
257  }
258  }
259  }
260  }
261  }
262  else
263  {
264  // For cell and vertex integrals use all markers (but not on ghost
265  // entities)
266  for (auto e = tagged_entities.begin(); e != entity_end; ++e)
267  {
268  const std::size_t i = std::distance(tagged_entities.cbegin(), e);
269  if (const auto it = id_to_integral.find(values[i]);
270  it != id_to_integral.end())
271  {
272  integrals[it->second].active_entities.push_back(*e);
273  }
274  }
275  }
276  }
277 
283  void set_default_domains(const mesh::Mesh& mesh)
284  {
285  const mesh::Topology& topology = mesh.topology();
286  const int tdim = topology.dim();
287  std::vector<struct Integral>& cell_integrals
288  = _integrals[static_cast<int>(IntegralType::cell)];
289 
290  // Cells. If there is a default integral, define it on all owned cells
291  if (cell_integrals.size() > 0 and cell_integrals[0].id == -1)
292  {
293  const int num_cells = topology.index_map(tdim)->size_local();
294  cell_integrals[0].active_entities.resize(num_cells);
295  std::iota(cell_integrals[0].active_entities.begin(),
296  cell_integrals[0].active_entities.end(), 0);
297  }
298 
299  // Exterior facets. If there is a default integral, define it only on
300  // owned surface facets.
301  std::vector<struct Integral>& exf_integrals
302  = _integrals[static_cast<int>(IntegralType::exterior_facet)];
303  if (exf_integrals.size() > 0 and exf_integrals[0].id == -1)
304  {
305  // If there is a default integral, define it only on surface facets
306  exf_integrals[0].active_entities.clear();
307 
308  // Get number of facets owned by this process
309  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
310  auto f_to_c = topology.connectivity(tdim - 1, tdim);
311  assert(topology.index_map(tdim - 1));
312  std::set<std::int32_t> fwd_shared_facets;
313 
314  // Only need to consider shared facets when there are no ghost cells
315  if (topology.index_map(tdim)->num_ghosts() == 0)
316  {
317  fwd_shared_facets.insert(
318  topology.index_map(tdim - 1)->shared_indices().begin(),
319  topology.index_map(tdim - 1)->shared_indices().end());
320  }
321 
322  const int num_facets = topology.index_map(tdim - 1)->size_local();
323  for (int f = 0; f < num_facets; ++f)
324  {
325  if (f_to_c->num_links(f) == 1
326  and fwd_shared_facets.find(f) == fwd_shared_facets.end())
327  exf_integrals[0].active_entities.push_back(f);
328  }
329  }
330 
331  // Interior facets. If there is a default integral, define it only on
332  // owned interior facets.
333  std::vector<struct FormIntegrals::Integral>& inf_integrals
334  = _integrals[static_cast<int>(IntegralType::interior_facet)];
335  if (!inf_integrals.empty() and inf_integrals[0].id == -1)
336  {
337  // If there is a default integral, define it only on interior facets
338  inf_integrals[0].active_entities.clear();
339 
340  // Get number of facets owned by this process
341  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
342  assert(topology.index_map(tdim - 1));
343  const int num_facets = topology.index_map(tdim - 1)->size_local();
344  auto f_to_c = topology.connectivity(tdim - 1, tdim);
345  inf_integrals[0].active_entities.reserve(num_facets);
346  for (int f = 0; f < num_facets; ++f)
347  {
348  if (f_to_c->num_links(f) == 2)
349  inf_integrals[0].active_entities.push_back(f);
350  }
351  }
352  }
353 
357  bool needs_permutation_data() const { return _needs_permutation_data; }
358 
359 private:
360  // Collect together the function, id, and indices of entities to
361  // integrate on
362  struct Integral
363  {
364  std::function<void(T*, const T*, const T*, const double*, const int*,
365  const std::uint8_t*, const std::uint32_t)>
366  tabulate;
367  int id;
368  std::vector<std::int32_t> active_entities;
369  };
370 
371  // Array of vectors of integrals, arranged by type (see Type enum, and
372  // struct Integral above)
373  std::array<std::vector<struct Integral>, 4> _integrals;
374 
375  // A bool indicating whether permutation data needs to be passed into
376  // these integrals
377  bool _needs_permutation_data;
378 };
379 } // namespace fem
380 } // namespace dolfinx
dolfinx::fem::FormIntegrals::num_integrals
int num_integrals(IntegralType type) const
Number of integrals of given type.
Definition: FormIntegrals.h:131
dolfinx::fem::FormIntegrals::set_domains
void set_domains(IntegralType type, const mesh::MeshTags< int > &marker)
Set the valid domains for the integrals of a given type from a MeshTags "marker". Note the MeshTags i...
Definition: FormIntegrals.h:168
dolfinx::fem::IntegralType
IntegralType
Type of integral.
Definition: FormIntegrals.h:27
dolfinx::fem::FormIntegrals::FormIntegrals
FormIntegrals(const std::map< IntegralType, std::vector< std::pair< int, std::function< void(T *, const T *, const T *, const double *, const int *, const std::uint8_t *, const std::uint32_t)>>>> &integrals, bool needs_permutation_data)
Construct object from (index, tabulate function) pairs for different integral types.
Definition: FormIntegrals.h:48
dolfinx::mesh::Mesh::topology_mutable
Topology & topology_mutable() const
Get mesh topology if one really needs the mutable version.
Definition: Mesh.cpp:145
dolfinx::fem::FormIntegrals::integral_domains
const std::vector< std::int32_t > & integral_domains(IntegralType type, int i) const
Get the list of active entities for the ith integral of type t. Note, these are not retrieved by ID,...
Definition: FormIntegrals.h:157
dolfinx::mesh::Topology::dim
int dim() const
Return topological dimension.
Definition: Topology.cpp:153
dolfinx::mesh::MeshTags::dim
int dim() const
Return topological dimension of tagged entities.
Definition: MeshTags.h:104
dolfinx::fem::FormIntegrals::types
std::set< IntegralType > types() const
Get types of integrals in the form.
Definition: FormIntegrals.h:116
dolfinx::fem::FormIntegrals::set_tabulate_tensor
void set_tabulate_tensor(IntegralType type, int i, const std::function< void(T *, const T *, const T *, const double *, const int *, const std::uint8_t *, const std::uint32_t)> &fn)
Definition: FormIntegrals.h:87
dolfinx::mesh::MeshTags::mesh
std::shared_ptr< const Mesh > mesh() const
Return mesh.
Definition: MeshTags.h:107
dolfinx::mesh::MeshTags
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: Form.h:36
dolfinx::fem::FormIntegrals::set_default_domains
void set_default_domains(const mesh::Mesh &mesh)
If there exists a default integral of any type, set the list of entities for those integrals from the...
Definition: FormIntegrals.h:283
dolfinx::mesh::Topology::index_map
std::shared_ptr< const common::IndexMap > index_map(int dim) const
Get the IndexMap that described the parallel distribution of the mesh entities.
Definition: Topology.cpp:162
dolfinx::fem::FormIntegrals::get_tabulate_tensor
const std::function< void(T *, const T *, const T *, const double *, const int *, const std::uint8_t *, const std::uint32_t)> & get_tabulate_tensor(IntegralType type, int i) const
Get the function for 'tabulate_tensor' for integral i of given type.
Definition: FormIntegrals.h:75
dolfinx::mesh::Mesh
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition: Mesh.h:47
dolfinx::mesh::Topology::connectivity
std::shared_ptr< const graph::AdjacencyList< std::int32_t > > connectivity(int d0, int d1) const
Return connectivity from entities of dimension d0 to entities of dimension d1.
Definition: Topology.cpp:255
dolfinx::mesh::Mesh::topology
Topology & topology()
Get mesh topology.
Definition: Mesh.cpp:141
dolfinx::mesh::MeshTags::indices
const std::vector< std::int32_t > & indices() const
Indices of tagged mesh entities (local-to-process). The indices are sorted.
Definition: MeshTags.h:98
dolfinx::mesh::MeshTags::values
const std::vector< T > & values() const
Values attached to mesh entities.
Definition: MeshTags.h:101
dolfinx::mesh::Topology::create_connectivity
void create_connectivity(int d0, int d1)
Create connectivity between given pair of dimensions, d0 -> d1.
Definition: Topology.cpp:194
dolfinx::fem::FormIntegrals::needs_permutation_data
bool needs_permutation_data() const
Get bool indicating whether permutation data needs to be passed into these integrals.
Definition: FormIntegrals.h:357
dolfinx::mesh::Topology
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition: Topology.h:58
dolfinx::fem::FormIntegrals
Integrals of a Form, including those defined over cells, interior and exterior facets,...
Definition: FormIntegrals.h:39
dolfinx::fem::FormIntegrals::integral_ids
std::vector< int > integral_ids(IntegralType type) const
Get the integer IDs of integrals of type t. The IDs correspond to the domains which the integrals are...
Definition: FormIntegrals.h:141