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:
43 
49  const std::function<void(T*, const T*, const T*, const double*, const int*,
50  const std::uint8_t*, const std::uint32_t)>&
51  get_tabulate_tensor(IntegralType type, int i) const
52  {
53  return _integrals.at(static_cast<int>(type)).at(i).tabulate;
54  }
55 
62  IntegralType type, int i,
63  std::function<void(T*, const T*, const T*, const double*, const int*,
64  const std::uint8_t*, const std::uint32_t)>
65  fn)
66  {
67  std::vector<struct FormIntegrals::Integral>& integrals
68  = _integrals.at(static_cast<int>(type));
69 
70  // Find insertion point
71  int pos = 0;
72  for (const auto& q : integrals)
73  {
74  if (q.id == i)
75  {
76  throw std::runtime_error("Integral with ID " + std::to_string(i)
77  + " already exists");
78  }
79  else if (q.id > i)
80  break;
81  ++pos;
82  }
83 
84  // Insert new Integral
85  integrals.insert(integrals.begin() + pos, {fn, i, {}});
86  }
87 
90  std::set<IntegralType> types() const
91  {
92  static const std::array types{IntegralType::cell,
93  IntegralType::exterior_facet,
94  IntegralType::interior_facet};
95  std::set<IntegralType> set;
96  for (auto type : types)
97  if (!_integrals.at(static_cast<int>(type)).empty())
98  set.insert(type);
99  return set;
100  }
101 
105  int num_integrals(IntegralType type) const
106  {
107  return _integrals.at(static_cast<int>(type)).size();
108  }
109 
115  std::vector<int> integral_ids(IntegralType type) const
116  {
117  std::vector<int> ids;
118  for (const auto& integral : _integrals[static_cast<int>(type)])
119  ids.push_back(integral.id);
120  return ids;
121  }
122 
131  const std::vector<std::int32_t>& integral_domains(IntegralType type,
132  int i) const
133  {
134  return _integrals.at(static_cast<int>(type)).at(i).active_entities;
135  }
136 
142  void set_domains(IntegralType type, const mesh::MeshTags<int>& marker)
143  {
144  std::vector<struct Integral>& integrals
145  = _integrals.at(static_cast<int>(type));
146  if (integrals.size() == 0)
147  return;
148 
149  std::shared_ptr<const mesh::Mesh> mesh = marker.mesh();
150  const mesh::Topology& topology = mesh->topology();
151  const int tdim = topology.dim();
152  int dim = tdim;
153  if (type == IntegralType::exterior_facet
154  or type == IntegralType::interior_facet)
155  {
156  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
157  dim = tdim - 1;
158  }
159  else if (type == IntegralType::vertex)
160  dim = 0;
161 
162  if (dim != marker.dim())
163  {
164  throw std::runtime_error("Invalid MeshTags dimension:"
165  + std::to_string(marker.dim()));
166  }
167 
168  // Create a reverse map
169  std::map<int, int> id_to_integral;
170  for (std::size_t i = 0; i < integrals.size(); ++i)
171  {
172  if (integrals[i].id != -1)
173  {
174  integrals[i].active_entities.clear();
175  id_to_integral.insert({integrals[i].id, i});
176  }
177  }
178 
179  // Get mesh tag data
180  const std::vector<int>& values = marker.values();
181  const std::vector<std::int32_t>& tagged_entities = marker.indices();
182  assert(topology.index_map(dim));
183  const auto entity_end
184  = std::lower_bound(tagged_entities.begin(), tagged_entities.end(),
185  topology.index_map(dim)->size_local());
186 
187  if (type == IntegralType::exterior_facet)
188  {
189  auto f_to_c = topology.connectivity(tdim - 1, tdim);
190  assert(f_to_c);
191  if (type == IntegralType::exterior_facet)
192  {
193  // Only need to consider shared facets when there are no ghost cells
194  assert(topology.index_map(tdim));
195  std::set<std::int32_t> fwd_shared;
196  if (topology.index_map(tdim)->num_ghosts() == 0)
197  {
198  fwd_shared.insert(
199  topology.index_map(tdim - 1)->shared_indices().begin(),
200  topology.index_map(tdim - 1)->shared_indices().end());
201  }
202 
203  for (auto f = tagged_entities.begin(); f != entity_end; ++f)
204  {
205  const std::size_t i = std::distance(tagged_entities.cbegin(), f);
206 
207  // All "owned" facets connected to one cell, that are not shared,
208  // should be external.
209  if (f_to_c->num_links(*f) == 1
210  and fwd_shared.find(*f) == fwd_shared.end())
211  {
212  if (auto it = id_to_integral.find(values[i]);
213  it != id_to_integral.end())
214  {
215  integrals[it->second].active_entities.push_back(*f);
216  }
217  }
218  }
219  }
220  else if (type == IntegralType::interior_facet)
221  {
222  for (auto f = tagged_entities.begin(); f != entity_end; ++f)
223  {
224  if (f_to_c->num_links(*f) == 2)
225  {
226  const std::size_t i = std::distance(tagged_entities.cbegin(), f);
227  if (const auto it = id_to_integral.find(values[i]);
228  it != id_to_integral.end())
229  {
230  integrals[it->second].active_entities.push_back(*f);
231  }
232  }
233  }
234  }
235  }
236  else
237  {
238  // For cell and vertex integrals use all markers (but not on ghost
239  // entities)
240  for (auto e = tagged_entities.begin(); e != entity_end; ++e)
241  {
242  const std::size_t i = std::distance(tagged_entities.cbegin(), e);
243  if (const auto it = id_to_integral.find(values[i]);
244  it != id_to_integral.end())
245  {
246  integrals[it->second].active_entities.push_back(*e);
247  }
248  }
249  }
250  }
251 
257  void set_default_domains(const mesh::Mesh& mesh)
258  {
259  const mesh::Topology& topology = mesh.topology();
260  const int tdim = topology.dim();
261  std::vector<struct Integral>& cell_integrals
262  = _integrals[static_cast<int>(IntegralType::cell)];
263 
264  // Cells. If there is a default integral, define it on all owned cells
265  if (cell_integrals.size() > 0 and cell_integrals[0].id == -1)
266  {
267  const int num_cells = topology.index_map(tdim)->size_local();
268  cell_integrals[0].active_entities.resize(num_cells);
269  std::iota(cell_integrals[0].active_entities.begin(),
270  cell_integrals[0].active_entities.end(), 0);
271  }
272 
273  // Exterior facets. If there is a default integral, define it only on
274  // owned surface facets.
275  std::vector<struct Integral>& exf_integrals
276  = _integrals[static_cast<int>(IntegralType::exterior_facet)];
277  if (exf_integrals.size() > 0 and exf_integrals[0].id == -1)
278  {
279  // If there is a default integral, define it only on surface facets
280  exf_integrals[0].active_entities.clear();
281 
282  // Get number of facets owned by this process
283  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
284  auto f_to_c = topology.connectivity(tdim - 1, tdim);
285  assert(topology.index_map(tdim - 1));
286  std::set<std::int32_t> fwd_shared_facets;
287 
288  // Only need to consider shared facets when there are no ghost cells
289  if (topology.index_map(tdim)->num_ghosts() == 0)
290  {
291  fwd_shared_facets.insert(
292  topology.index_map(tdim - 1)->shared_indices().begin(),
293  topology.index_map(tdim - 1)->shared_indices().end());
294  }
295 
296  const int num_facets = topology.index_map(tdim - 1)->size_local();
297  for (int f = 0; f < num_facets; ++f)
298  {
299  if (f_to_c->num_links(f) == 1
300  and fwd_shared_facets.find(f) == fwd_shared_facets.end())
301  exf_integrals[0].active_entities.push_back(f);
302  }
303  }
304 
305  // Interior facets. If there is a default integral, define it only on
306  // owned interior facets.
307  std::vector<struct FormIntegrals::Integral>& inf_integrals
308  = _integrals[static_cast<int>(IntegralType::interior_facet)];
309  if (!inf_integrals.empty() and inf_integrals[0].id == -1)
310  {
311  // If there is a default integral, define it only on interior facets
312  inf_integrals[0].active_entities.clear();
313 
314  // Get number of facets owned by this process
315  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
316  assert(topology.index_map(tdim - 1));
317  const int num_facets = topology.index_map(tdim - 1)->size_local();
318  auto f_to_c = topology.connectivity(tdim - 1, tdim);
319  inf_integrals[0].active_entities.reserve(num_facets);
320  for (int f = 0; f < num_facets; ++f)
321  {
322  if (f_to_c->num_links(f) == 2)
323  inf_integrals[0].active_entities.push_back(f);
324  }
325  }
326  }
327 
328 private:
329  // Collect together the function, id, and indices of entities to
330  // integrate on
331  struct Integral
332  {
333  std::function<void(T*, const T*, const T*, const double*, const int*,
334  const std::uint8_t*, const std::uint32_t)>
335  tabulate;
336  int id;
337  std::vector<std::int32_t> active_entities;
338  };
339 
340  // Array of vectors of integrals, arranged by type (see Type enum, and
341  // struct Integral above)
342  std::array<std::vector<struct Integral>, 4> _integrals;
343 };
344 } // namespace fem
345 } // namespace dolfinx
dolfinx::fem::FormIntegrals::num_integrals
int num_integrals(IntegralType type) const
Number of integrals of given type.
Definition: FormIntegrals.h:105
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:142
dolfinx::fem::IntegralType
IntegralType
Type of integral.
Definition: FormIntegrals.h:27
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:131
dolfinx::mesh::Topology::dim
int dim() const
Return topological dimension.
Definition: Topology.cpp:152
dolfinx::mesh::MeshTags::dim
int dim() const
Return topological dimension of tagged entities.
Definition: MeshTags.h:104
dolfinx::fem::FormIntegrals::FormIntegrals
FormIntegrals()
Construct empty object.
Definition: FormIntegrals.h:42
dolfinx::fem::FormIntegrals::types
std::set< IntegralType > types() const
Get types of integrals in the form.
Definition: FormIntegrals.h:90
dolfinx::mesh::MeshTags::mesh
std::shared_ptr< const Mesh > mesh() const
Return mesh.
Definition: MeshTags.h:107
dolfinx::fem::FormIntegrals::set_tabulate_tensor
void set_tabulate_tensor(IntegralType type, int i, std::function< void(T *, const T *, const T *, const double *, const int *, const std::uint8_t *, const std::uint32_t)> fn)
Set the function for 'tabulate_tensor' for integral i of given type.
Definition: FormIntegrals.h:61
dolfinx::mesh::MeshTags
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: Form.h:37
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:257
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:161
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:51
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:254
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:193
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:115