DOLFIN-X
DOLFIN-X C++ interface
Form.h
1 // Copyright (C) 2019-2020 Garth N. Wells and 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 <dolfinx/function/FunctionSpace.h>
10 #include <dolfinx/mesh/Mesh.h>
11 #include <dolfinx/mesh/MeshTags.h>
12 #include <functional>
13 #include <memory>
14 #include <string>
15 #include <vector>
16 
17 namespace dolfinx
18 {
19 
20 namespace function
21 {
22 template <typename T>
23 class Constant;
24 template <typename T>
25 class Function;
26 } // namespace function
27 
28 namespace fem
29 {
30 
32 enum class IntegralType : std::int8_t
33 {
34  cell = 0,
35  exterior_facet = 1,
36  interior_facet = 2,
37  vertex = 3
38 };
39 
63 
64 template <typename T>
65 class Form
66 {
67 public:
81  Form(const std::vector<std::shared_ptr<const function::FunctionSpace>>&
83  const std::map<
85  std::pair<
86  std::vector<std::pair<
87  int, std::function<void(
88  T*, const T*, const T*, const double*, const int*,
89  const std::uint8_t*, const std::uint32_t)>>>,
90  const mesh::MeshTags<int>*>>& integrals,
91  const std::vector<std::shared_ptr<const function::Function<T>>>&
93  const std::vector<std::shared_ptr<const function::Constant<T>>>&
94  constants,
96  const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
97  : _function_spaces(function_spaces), _coefficients(coefficients),
98  _constants(constants), _mesh(mesh),
99  _needs_permutation_data(needs_permutation_data)
100  {
101  // Extract _mesh from function::FunctionSpace, and check they are the same
102  if (!_mesh and !function_spaces.empty())
103  _mesh = function_spaces[0]->mesh();
104  for (const auto& V : function_spaces)
105  if (_mesh != V->mesh())
106  throw std::runtime_error("Incompatible mesh");
107  if (!_mesh)
108  throw std::runtime_error("No mesh could be associated with the Form.");
109 
110  // Store kernels, looping over integrals by domain type (dimension)
111  for (auto& integral_type : integrals)
112  {
113  // Add key to map
114  const IntegralType type = integral_type.first;
115  auto it = _integrals.emplace(
116  type, std::map<int, std::pair<kern, std::vector<std::int32_t>>>());
117 
118  // Loop over integrals kernels
119  for (auto& integral : integral_type.second.first)
120  it.first->second.insert({integral.first, {integral.second, {}}});
121 
122  // FIXME: do this neatly via a static function
123  // Set domains for integral type
124  if (integral_type.second.second)
125  {
126  assert(_mesh == integral_type.second.second->mesh());
127  set_domains(type, *integral_type.second.second);
128  }
129  }
130 
131  // FIXME: do this neatly via a static function
132  // Set markers for default integrals
133  set_default_domains(*_mesh);
134  }
135 
137  Form(const Form& form) = delete;
138 
140  Form(Form&& form) = default;
141 
143  virtual ~Form() = default;
144 
148  int rank() const { return _function_spaces.size(); }
149 
152  std::shared_ptr<const mesh::Mesh> mesh() const { return _mesh; }
153 
156  const std::vector<std::shared_ptr<const function::FunctionSpace>>&
158  {
159  return _function_spaces;
160  }
161 
167  const std::function<void(T*, const T*, const T*, const double*, const int*,
168  const std::uint8_t*, const std::uint32_t)>&
169  kernel(IntegralType type, int i) const
170  {
171  auto it0 = _integrals.find(type);
172  if (it0 == _integrals.end())
173  throw std::runtime_error("No kernels for requested type.");
174  auto it1 = it0->second.find(i);
175  if (it1 == it0->second.end())
176  throw std::runtime_error("No kernel for requested domain index.");
177 
178  return it1->second.first;
179  }
180 
183  std::set<IntegralType> integral_types() const
184  {
185  std::set<IntegralType> set;
186  for (auto& type : _integrals)
187  set.insert(type.first);
188  return set;
189  }
190 
194  int num_integrals(IntegralType type) const
195  {
196  if (auto it = _integrals.find(type); it == _integrals.end())
197  return 0;
198  else
199  return it->second.size();
200  }
201 
208  std::vector<int> integral_ids(IntegralType type) const
209  {
210  std::vector<int> ids;
211  if (auto it = _integrals.find(type); it != _integrals.end())
212  {
213  for (auto& kernel : it->second)
214  ids.push_back(kernel.first);
215  }
216  return ids;
217  }
218 
225  const std::vector<std::int32_t>& domains(IntegralType type, int i) const
226  {
227  auto it0 = _integrals.find(type);
228  if (it0 == _integrals.end())
229  throw std::runtime_error("No mesh entities for requested type.");
230  auto it1 = it0->second.find(i);
231  if (it1 == it0->second.end())
232  throw std::runtime_error("No mesh entities for requested domain index.");
233  return it1->second.second;
234  }
235 
237  const std::vector<std::shared_ptr<const function::Function<T>>>
238  coefficients() const
239  {
240  return _coefficients;
241  }
242 
246  bool needs_permutation_data() const { return _needs_permutation_data; }
247 
251  std::vector<int> coefficient_offsets() const
252  {
253  std::vector<int> n{0};
254  for (const auto& c : _coefficients)
255  {
256  if (!c)
257  throw std::runtime_error("Not all form coefficients have been set.");
258  n.push_back(n.back() + c->function_space()->element()->space_dimension());
259  }
260  return n;
261  }
262 
264  const std::vector<std::shared_ptr<const function::Constant<T>>>&
265  constants() const
266  {
267  return _constants;
268  }
269 
271  using scalar_type = T;
272 
273 private:
279  void set_domains(IntegralType type, const mesh::MeshTags<int>& marker)
280  {
281  auto it0 = _integrals.find(type);
282  assert(it0 != _integrals.end());
283 
284  std::shared_ptr<const mesh::Mesh> mesh = marker.mesh();
285  const mesh::Topology& topology = mesh->topology();
286  const int tdim = topology.dim();
287  int dim = tdim;
288  if (type == IntegralType::exterior_facet
289  or type == IntegralType::interior_facet)
290  {
291  dim = tdim - 1;
292  mesh->topology_mutable().create_connectivity(dim, tdim);
293  }
294  else if (type == IntegralType::vertex)
295  dim = 0;
296 
297  if (dim != marker.dim())
298  {
299  throw std::runtime_error("Invalid MeshTags dimension:"
300  + std::to_string(marker.dim()));
301  }
302 
303  // Get all integrals for considered entity type
304  std::map<int, std::pair<kern, std::vector<std::int32_t>>>& integrals
305  = it0->second;
306 
307  // Get mesh tag data
308  const std::vector<int>& values = marker.values();
309  const std::vector<std::int32_t>& tagged_entities = marker.indices();
310  assert(topology.index_map(dim));
311  const auto entity_end
312  = std::lower_bound(tagged_entities.begin(), tagged_entities.end(),
313  topology.index_map(dim)->size_local());
314 
315  if (dim == tdim - 1)
316  {
317  auto f_to_c = topology.connectivity(tdim - 1, tdim);
318  assert(f_to_c);
319  if (type == IntegralType::exterior_facet)
320  {
321  // Only need to consider shared facets when there are no ghost
322  // cells
323  assert(topology.index_map(tdim));
324  std::set<std::int32_t> fwd_shared;
325  if (topology.index_map(tdim)->num_ghosts() == 0)
326  {
327  fwd_shared.insert(
328  topology.index_map(tdim - 1)->shared_indices().begin(),
329  topology.index_map(tdim - 1)->shared_indices().end());
330  }
331 
332  for (auto f = tagged_entities.begin(); f != entity_end; ++f)
333  {
334  // All "owned" facets connected to one cell, that are not
335  // shared, should be external
336  if (f_to_c->num_links(*f) == 1
337  and fwd_shared.find(*f) == fwd_shared.end())
338  {
339  const std::size_t i = std::distance(tagged_entities.cbegin(), f);
340  if (auto it = integrals.find(values[i]); it != integrals.end())
341  it->second.second.push_back(*f);
342  }
343  }
344  }
345  else if (type == IntegralType::interior_facet)
346  {
347  for (auto f = tagged_entities.begin(); f != entity_end; ++f)
348  {
349  if (f_to_c->num_links(*f) == 2)
350  {
351  const std::size_t i = std::distance(tagged_entities.cbegin(), f);
352  if (auto it = integrals.find(values[i]); it != integrals.end())
353  it->second.second.push_back(*f);
354  }
355  }
356  }
357  }
358  else
359  {
360  // For cell and vertex integrals use all markers (but not on ghost
361  // entities)
362  for (auto e = tagged_entities.begin(); e != entity_end; ++e)
363  {
364  const std::size_t i = std::distance(tagged_entities.cbegin(), e);
365  if (auto it = integrals.find(values[i]); it != integrals.end())
366  it->second.second.push_back(*e);
367  }
368  }
369  }
370 
376  void set_default_domains(const mesh::Mesh& mesh)
377  {
378  const mesh::Topology& topology = mesh.topology();
379  const int tdim = topology.dim();
380 
381  // Cells. If there is a default integral, define it on all owned cells
382  if (auto kernels = _integrals.find(IntegralType::cell);
383  kernels != _integrals.end())
384  {
385  if (auto it = kernels->second.find(-1); it != kernels->second.end())
386  {
387  std::vector<std::int32_t>& active_entities = it->second.second;
388  const int num_cells = topology.index_map(tdim)->size_local();
389  active_entities.resize(num_cells);
390  std::iota(active_entities.begin(), active_entities.end(), 0);
391  }
392  }
393 
394  // Exterior facets. If there is a default integral, define it only
395  // on owned surface facets.
396  if (auto kernels = _integrals.find(IntegralType::exterior_facet);
397  kernels != _integrals.end())
398  {
399  if (auto it = kernels->second.find(-1); it != kernels->second.end())
400  {
401  std::vector<std::int32_t>& active_entities = it->second.second;
402  active_entities.clear();
403 
404  // Get number of facets owned by this process
405  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
406  auto f_to_c = topology.connectivity(tdim - 1, tdim);
407  assert(topology.index_map(tdim - 1));
408  std::set<std::int32_t> fwd_shared_facets;
409 
410  // Only need to consider shared facets when there are no ghost cells
411  if (topology.index_map(tdim)->num_ghosts() == 0)
412  {
413  fwd_shared_facets.insert(
414  topology.index_map(tdim - 1)->shared_indices().begin(),
415  topology.index_map(tdim - 1)->shared_indices().end());
416  }
417 
418  const int num_facets = topology.index_map(tdim - 1)->size_local();
419  for (int f = 0; f < num_facets; ++f)
420  {
421  if (f_to_c->num_links(f) == 1
422  and fwd_shared_facets.find(f) == fwd_shared_facets.end())
423  {
424  active_entities.push_back(f);
425  }
426  }
427  }
428  }
429 
430  // Interior facets. If there is a default integral, define it only on
431  // owned interior facets.
432  if (auto kernels = _integrals.find(IntegralType::interior_facet);
433  kernels != _integrals.end())
434  {
435  if (auto it = kernels->second.find(-1); it != kernels->second.end())
436  {
437  std::vector<std::int32_t>& active_entities = it->second.second;
438 
439  // Get number of facets owned by this process
440  mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
441  assert(topology.index_map(tdim - 1));
442  const int num_facets = topology.index_map(tdim - 1)->size_local();
443  auto f_to_c = topology.connectivity(tdim - 1, tdim);
444  active_entities.clear();
445  active_entities.reserve(num_facets);
446  for (int f = 0; f < num_facets; ++f)
447  {
448  if (f_to_c->num_links(f) == 2)
449  active_entities.push_back(f);
450  }
451  }
452  }
453  }
454 
455  // Function spaces (one for each argument)
456  std::vector<std::shared_ptr<const function::FunctionSpace>> _function_spaces;
457 
458  // Form coefficients
459  std::vector<std::shared_ptr<const function::Function<T>>> _coefficients;
460 
461  // Constants associated with the Form
462  std::vector<std::shared_ptr<const function::Constant<T>>> _constants;
463 
464  // The mesh
465  std::shared_ptr<const mesh::Mesh> _mesh;
466 
467  using kern
468  = std::function<void(T*, const T*, const T*, const double*, const int*,
469  const std::uint8_t*, const std::uint32_t)>;
470  std::map<IntegralType,
471  std::map<int, std::pair<kern, std::vector<std::int32_t>>>>
472  _integrals;
473 
474  // True if permutation data needs to be passed into these integrals
475  bool _needs_permutation_data;
476 
477 }; // namespace fem
478 
479 } // namespace fem
480 } // namespace dolfinx
dolfinx::fem::Form::coefficient_offsets
std::vector< int > coefficient_offsets() const
Offset for each coefficient expansion array on a cell. Used to pack data for multiple coefficients in...
Definition: Form.h:251
dolfinx::fem::IntegralType
IntegralType
Type of integral.
Definition: Form.h:33
dolfinx::fem::Form::constants
const std::vector< std::shared_ptr< const function::Constant< T > > > & constants() const
Access constants.
Definition: Form.h:265
dolfinx::fem::Form::integral_types
std::set< IntegralType > integral_types() const
Get types of integrals in the form.
Definition: Form.h:183
dolfinx::fem::Form::rank
int rank() const
Rank of the form (bilinear form = 2, linear form = 1, functional = 0, etc)
Definition: Form.h:148
dolfinx::fem::Form::num_integrals
int num_integrals(IntegralType type) const
Number of integrals of given type.
Definition: Form.h:194
dolfinx::fem::Form::function_spaces
const std::vector< std::shared_ptr< const function::FunctionSpace > > & function_spaces() const
Return function spaces for all arguments.
Definition: Form.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::Form::~Form
virtual ~Form()=default
Destructor.
dolfinx::fem::Form
Class for variational forms.
Definition: Form.h:66
dolfinx::mesh::MeshTags::mesh
std::shared_ptr< const Mesh > mesh() const
Return mesh.
Definition: MeshTags.h:107
dolfinx::fem::Form::Form
Form(Form &&form)=default
Move constructor.
dolfinx::mesh::MeshTags
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: MeshTags.h:37
dolfinx::fem::Form::Form
Form(const Form &form)=delete
Copy constructor.
dolfinx::function::Constant
A constant value which can be attached to a Form. Constants may be scalar (rank 0),...
Definition: Constant.h:19
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::Form::Form
Form(const std::vector< std::shared_ptr< const function::FunctionSpace >> &function_spaces, const std::map< IntegralType, std::pair< 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)>>>, const mesh::MeshTags< int > * >> &integrals, const std::vector< std::shared_ptr< const function::Function< T >>> &coefficients, const std::vector< std::shared_ptr< const function::Constant< T >>> &constants, bool needs_permutation_data, const std::shared_ptr< const mesh::Mesh > &mesh=nullptr)
Create form.
Definition: Form.h:81
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::fem::Form::kernel
const std::function< void(T *, const T *, const T *, const double *, const int *, const std::uint8_t *, const std::uint32_t)> & kernel(IntegralType type, int i) const
Get the function for 'kernel' for integral i of given type.
Definition: Form.h:169
dolfinx::fem::Form::integral_ids
std::vector< int > integral_ids(IntegralType type) const
Get the IDs for integrals (kernels) for given integral type. The IDs correspond to the domain IDs whi...
Definition: Form.h:208
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::fem::Form::mesh
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh for the form.
Definition: Form.h:152
dolfinx::mesh::MeshTags::values
const std::vector< T > & values() const
Values attached to mesh entities.
Definition: MeshTags.h:101
dolfinx::fem::Form::scalar_type
T scalar_type
Scalar type (T).
Definition: Form.h:271
dolfinx::fem::Form::needs_permutation_data
bool needs_permutation_data() const
Get bool indicating whether permutation data needs to be passed into these integrals.
Definition: Form.h:246
dolfinx::fem::Form::domains
const std::vector< std::int32_t > & domains(IntegralType type, int i) const
Get the list of mesh entity indices for the ith integral (kernel) for the given domain type,...
Definition: Form.h:225
dolfinx::function::Function
This class represents a function in a finite element function space , given by.
Definition: Function.h:42
dolfinx::mesh::Topology
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition: Topology.h:58
dolfinx::fem::Form::coefficients
const std::vector< std::shared_ptr< const function::Function< T > > > coefficients() const
Access coefficients.
Definition: Form.h:238