DOLFIN-X
DOLFIN-X C++ interface
utils.h
1 // Copyright (C) 2013-2020 Johan Hake, Jan Blechta and 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 "CoordinateElement.h"
10 #include "DofMap.h"
11 #include "ElementDofLayout.h"
12 #include <dolfinx/common/types.h>
13 #include <dolfinx/fem/Form.h>
14 #include <dolfinx/fem/Function.h>
15 #include <dolfinx/la/SparsityPattern.h>
16 #include <dolfinx/mesh/cell_types.h>
17 #include <memory>
18 #include <set>
19 #include <string>
20 #include <ufc.h>
21 #include <utility>
22 #include <vector>
23 
24 namespace dolfinx
25 {
26 namespace common
27 {
28 class IndexMap;
29 }
30 
31 namespace mesh
32 {
33 class Mesh;
34 class Topology;
35 } // namespace mesh
36 
37 namespace fem
38 {
39 
40 template <typename T>
41 class Constant;
42 template <typename T>
43 class Form;
44 template <typename T>
45 class Function;
46 class FunctionSpace;
47 
55 template <typename T>
56 std::vector<
57  std::vector<std::array<std::shared_ptr<const fem::FunctionSpace>, 2>>>
58 extract_function_spaces(const std::vector<std::vector<const fem::Form<T>*>>& a)
59 {
60  std::vector<
61  std::vector<std::array<std::shared_ptr<const fem::FunctionSpace>, 2>>>
62  spaces(
63  a.size(),
64  std::vector<std::array<std::shared_ptr<const fem::FunctionSpace>, 2>>(
65  a[0].size()));
66  for (std::size_t i = 0; i < a.size(); ++i)
67  {
68  for (std::size_t j = 0; j < a[i].size(); ++j)
69  {
70  if (const fem::Form<T>* form = a[i][j]; form)
71  spaces[i][j] = {form->function_spaces()[0], form->function_spaces()[1]};
72  }
73  }
74  return spaces;
75 }
76 
82 template <typename T>
84 {
85  if (a.rank() != 2)
86  {
87  throw std::runtime_error(
88  "Cannot create sparsity pattern. Form is not a bilinear form");
89  }
90 
91  // Get dof maps and mesh
92  std::array<const std::reference_wrapper<const fem::DofMap>, 2> dofmaps{
93  *a.function_spaces().at(0)->dofmap(),
94  *a.function_spaces().at(1)->dofmap()};
95  std::shared_ptr mesh = a.mesh();
96  assert(mesh);
97 
98  const std::set<IntegralType> types = a.integral_types();
99  if (types.find(IntegralType::interior_facet) != types.end()
100  or types.find(IntegralType::exterior_facet) != types.end())
101  {
102  // FIXME: cleanup these calls? Some of the happen internally again.
103  const int tdim = mesh->topology().dim();
104  mesh->topology_mutable().create_entities(tdim - 1);
105  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
106  }
107 
108  return create_sparsity_pattern(mesh->topology(), dofmaps, types);
109 }
110 
115  const mesh::Topology& topology,
116  const std::array<const std::reference_wrapper<const fem::DofMap>, 2>&
117  dofmaps,
118  const std::set<IntegralType>& integrals);
119 
121 ElementDofLayout create_element_dof_layout(const ufc_dofmap& dofmap,
122  const mesh::CellType cell_type,
123  const std::vector<int>& parent_map
124  = {});
125 
130 DofMap create_dofmap(MPI_Comm comm, const ufc_dofmap& dofmap,
131  mesh::Topology& topology);
132 
136 std::vector<std::string> get_coefficient_names(const ufc_form& ufc_form);
137 
141 std::vector<std::string> get_constant_names(const ufc_form& ufc_form);
142 
150 template <typename T>
152  const ufc_form& ufc_form,
153  const std::vector<std::shared_ptr<const fem::FunctionSpace>>& spaces,
154  const std::vector<std::shared_ptr<const fem::Function<T>>>& coefficients,
155  const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants,
156  const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
157  const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
158 {
159  if (ufc_form.rank != (int)spaces.size())
160  throw std::runtime_error("Wrong number of argument spaces for Form.");
161  if (ufc_form.num_coefficients != (int)coefficients.size())
162  {
163  throw std::runtime_error(
164  "Mismatch between number of expected and provided Form coefficients.");
165  }
166  if (ufc_form.num_constants != (int)constants.size())
167  {
168  throw std::runtime_error(
169  "Mismatch between number of expected and provided Form constants.");
170  }
171 
172  // Check argument function spaces
173 #ifdef DEBUG
174  for (std::size_t i = 0; i < spaces.size(); ++i)
175  {
176  assert(spaces[i]->element());
177  std::unique_ptr<ufc_finite_element, decltype(free)*> ufc_element(
178  ufc_form.create_finite_element(i), free);
179  assert(ufc_element);
180  if (std::string(ufc_element->signature)
181  != spaces[i]->element()->signature())
182  {
183  throw std::runtime_error(
184  "Cannot create form. Wrong type of function space for argument.");
185  }
186  }
187 #endif
188 
189  // Get list of integral IDs, and load tabulate tensor into memory for
190  // each
191  using kern
192  = std::function<void(T*, const T*, const T*, const double*, const int*,
193  const std::uint8_t*, const std::uint32_t)>;
194  std::map<IntegralType, std::pair<std::vector<std::pair<int, kern>>,
195  const mesh::MeshTags<int>*>>
196  integral_data;
197 
198  bool needs_permutation_data = false;
199 
200  // Attach cell kernels
201  std::vector<int> cell_integral_ids(ufc_form.num_cell_integrals);
202  ufc_form.get_cell_integral_ids(cell_integral_ids.data());
203  for (int id : cell_integral_ids)
204  {
205  ufc_integral* integral = ufc_form.create_cell_integral(id);
206  assert(integral);
207  if (integral->needs_permutation_data)
208  needs_permutation_data = true;
209  integral_data[IntegralType::cell].first.emplace_back(
210  id, integral->tabulate_tensor);
211  std::free(integral);
212  }
213 
214  // Attach cell subdomain data
215  if (auto it = subdomains.find(IntegralType::cell);
216  it != subdomains.end() and !cell_integral_ids.empty())
217  {
218  integral_data[IntegralType::cell].second = it->second;
219  }
220 
221  // FIXME: Can facets be handled better?
222 
223  // Create facets, if required
224  if (ufc_form.num_exterior_facet_integrals > 0
225  or ufc_form.num_interior_facet_integrals > 0)
226  {
227  if (!spaces.empty())
228  {
229  auto mesh = spaces[0]->mesh();
230  const int tdim = mesh->topology().dim();
231  spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
232  }
233  }
234 
235  // Attach exterior facet kernels
236  std::vector<int> exterior_facet_integral_ids(
237  ufc_form.num_exterior_facet_integrals);
238  ufc_form.get_exterior_facet_integral_ids(exterior_facet_integral_ids.data());
239  for (int id : exterior_facet_integral_ids)
240  {
241  ufc_integral* integral = ufc_form.create_exterior_facet_integral(id);
242  assert(integral);
243  if (integral->needs_permutation_data)
244  needs_permutation_data = true;
245  integral_data[IntegralType::exterior_facet].first.emplace_back(
246  id, integral->tabulate_tensor);
247  std::free(integral);
248  }
249 
250  // Attach exterior facet subdomain data
251  if (auto it = subdomains.find(IntegralType::exterior_facet);
252  it != subdomains.end() and !exterior_facet_integral_ids.empty())
253  {
254  integral_data[IntegralType::exterior_facet].second = it->second;
255  }
256 
257  // Attach interior facet kernels
258  std::vector<int> interior_facet_integral_ids(
259  ufc_form.num_interior_facet_integrals);
260  ufc_form.get_interior_facet_integral_ids(interior_facet_integral_ids.data());
261  for (int id : interior_facet_integral_ids)
262  {
263  ufc_integral* integral = ufc_form.create_interior_facet_integral(id);
264  assert(integral);
265  if (integral->needs_permutation_data)
266  needs_permutation_data = true;
267  integral_data[IntegralType::interior_facet].first.emplace_back(
268  id, integral->tabulate_tensor);
269  std::free(integral);
270  }
271 
272  // Attach interior facet subdomain data
273  if (auto it = subdomains.find(IntegralType::interior_facet);
274  it != subdomains.end() and !interior_facet_integral_ids.empty())
275  {
276  integral_data[IntegralType::interior_facet].second = it->second;
277  }
278 
279  // Vertex integrals: not currently working
280  std::vector<int> vertex_integral_ids(ufc_form.num_vertex_integrals);
281  ufc_form.get_vertex_integral_ids(vertex_integral_ids.data());
282  if (!vertex_integral_ids.empty())
283  {
284  throw std::runtime_error(
285  "Vertex integrals not supported. Under development.");
286  }
287 
288  return fem::Form(spaces, integral_data, coefficients, constants,
289  needs_permutation_data, mesh);
290 }
291 
301 template <typename T>
303  const ufc_form& ufc_form,
304  const std::vector<std::shared_ptr<const fem::FunctionSpace>>& spaces,
305  const std::map<std::string, std::shared_ptr<const fem::Function<T>>>&
306  coefficients,
307  const std::map<std::string, std::shared_ptr<const fem::Constant<T>>>&
308  constants,
309  const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
310  const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
311 {
312  // Place coefficients in appropriate order
313  std::vector<std::shared_ptr<const fem::Function<T>>> coeff_map;
314  for (const std::string& name : get_coefficient_names(ufc_form))
315  {
316  if (auto it = coefficients.find(name); it != coefficients.end())
317  coeff_map.push_back(it->second);
318  else
319  {
320  throw std::runtime_error("Form coefficient \"" + name
321  + "\" not provided.");
322  }
323  }
324 
325  // Place constants in appropriate order
326  std::vector<std::shared_ptr<const fem::Constant<T>>> const_map;
327  for (const std::string& name : get_constant_names(ufc_form))
328  {
329  if (auto it = constants.find(name); it != constants.end())
330  const_map.push_back(it->second);
331  else
332  {
333  throw std::runtime_error("Form constant \"" + name + "\" not provided.");
334  }
335  }
336 
337  return create_form(ufc_form, spaces, coeff_map, const_map, subdomains, mesh);
338 }
339 
351 template <typename T>
352 std::shared_ptr<Form<T>> create_form(
353  ufc_form* (*fptr)(),
354  const std::vector<std::shared_ptr<const fem::FunctionSpace>>& spaces,
355  const std::map<std::string, std::shared_ptr<const fem::Function<T>>>&
356  coefficients,
357  const std::map<std::string, std::shared_ptr<const fem::Constant<T>>>&
358  constants,
359  const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
360  const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
361 {
362  ufc_form* form = fptr();
363  auto L = std::make_shared<fem::Form<T>>(dolfinx::fem::create_form<T>(
364  *form, spaces, coefficients, constants, subdomains, mesh));
365  std::free(form);
366  return L;
367 }
368 
373 create_coordinate_map(const ufc_coordinate_mapping& ufc_cmap);
374 
379 fem::CoordinateElement create_coordinate_map(ufc_coordinate_mapping* (*fptr)());
380 
390 std::shared_ptr<fem::FunctionSpace>
391 create_functionspace(ufc_function_space* (*fptr)(const char*),
392  const std::string function_name,
393  std::shared_ptr<mesh::Mesh> mesh);
394 
395 // NOTE: This is subject to change
397 template <typename U>
398 Eigen::Array<typename U::scalar_type, Eigen::Dynamic, Eigen::Dynamic,
399  Eigen::RowMajor>
400 pack_coefficients(const U& u)
401 {
402  using T = typename U::scalar_type;
403 
404  // Get form coefficient offsets and dofmaps
405  const std::vector<std::shared_ptr<const fem::Function<T>>> coefficients
406  = u.coefficients();
407  const std::vector<int> offsets = u.coefficient_offsets();
408  std::vector<const fem::DofMap*> dofmaps(coefficients.size());
409  std::vector<int> bs(coefficients.size());
410  std::vector<std::reference_wrapper<const std::vector<T>>> v;
411  for (std::size_t i = 0; i < coefficients.size(); ++i)
412  {
413  dofmaps[i] = coefficients[i]->function_space()->dofmap().get();
414  bs[i] = dofmaps[i]->bs();
415  v.push_back(coefficients[i]->x()->array());
416  }
417 
418  // Get mesh
419  std::shared_ptr<const mesh::Mesh> mesh = u.mesh();
420  assert(mesh);
421  const int tdim = mesh->topology().dim();
422  const std::int32_t num_cells
423  = mesh->topology().index_map(tdim)->size_local()
424  + mesh->topology().index_map(tdim)->num_ghosts();
425 
426  // Copy data into coefficient array
427  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
428  num_cells, offsets.back());
429  if (coefficients.size() > 0)
430  {
431  for (int cell = 0; cell < num_cells; ++cell)
432  {
433  for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
434  {
435  tcb::span<const std::int32_t> dofs = dofmaps[coeff]->cell_dofs(cell);
436  const std::vector<T>& _v = v[coeff];
437  for (std::size_t i = 0; i < dofs.size(); ++i)
438  {
439  for (int k = 0; k < bs[coeff]; ++k)
440  {
441  c(cell, bs[coeff] * i + k + offsets[coeff])
442  = _v[bs[coeff] * dofs[i] + k];
443  }
444  }
445  }
446  }
447  }
448 
449  return c;
450 }
451 
452 // NOTE: This is subject to change
454 template <typename U>
455 std::vector<typename U::scalar_type> pack_constants(const U& u)
456 {
457  using T = typename U::scalar_type;
458  const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants
459  = u.constants();
460 
461  // Calculate size of array needed to store packed constants
462  std::int32_t size
463  = std::accumulate(constants.begin(), constants.end(), 0,
464  [](std::int32_t sum, const auto& constant) {
465  return sum + constant->value.size();
466  });
467 
468  // Pack constants
469  std::vector<T> constant_values(size);
470  std::int32_t offset = 0;
471  for (const auto& constant : constants)
472  {
473  const std::vector<T>& value = constant->value;
474  for (std::size_t i = 0; i < value.size(); ++i)
475  constant_values[offset + i] = value[i];
476  offset += value.size();
477  }
478 
479  return constant_values;
480 }
481 
482 } // namespace fem
483 } // namespace dolfinx
A constant value which can be attached to a Form. Constants may be scalar (rank 0),...
Definition: Constant.h:19
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:26
Class for variational forms.
Definition: Form.h:60
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh for the form.
Definition: Form.h:144
int rank() const
Rank of the form (bilinear form = 2, linear form = 1, functional = 0, etc)
Definition: Form.h:140
std::set< IntegralType > integral_types() const
Get types of integrals in the form.
Definition: Form.h:175
const std::vector< std::shared_ptr< const fem::FunctionSpace > > & function_spaces() const
Return function spaces for all arguments.
Definition: Form.h:149
This class represents a function in a finite element function space , given by.
Definition: Function.h:43
This class provides a sparsity pattern data structure that can be used to initialize sparse matrices.
Definition: SparsityPattern.h:37
A MeshTags are used to associate mesh entities with values. The entity index (local to process) ident...
Definition: MeshTags.h:38
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition: Topology.h:57
std::vector< typename U::scalar_type > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:455
Form< T > create_form(const ufc_form &ufc_form, const std::vector< std::shared_ptr< const fem::FunctionSpace >> &spaces, const std::vector< std::shared_ptr< const fem::Function< T >>> &coefficients, const std::vector< std::shared_ptr< const fem::Constant< T >>> &constants, const std::map< IntegralType, const mesh::MeshTags< int > * > &subdomains, const std::shared_ptr< const mesh::Mesh > &mesh=nullptr)
Create a Form from UFC input.
Definition: utils.h:151
ElementDofLayout create_element_dof_layout(const ufc_dofmap &dofmap, const mesh::CellType cell_type, const std::vector< int > &parent_map={})
Create an ElementDofLayout from a ufc_dofmap.
Definition: utils.cpp:76
std::shared_ptr< fem::FunctionSpace > create_functionspace(ufc_function_space *(*fptr)(const char *), const std::string function_name, std::shared_ptr< mesh::Mesh > mesh)
Create FunctionSpace from UFC.
Definition: utils.cpp:245
DofMap create_dofmap(MPI_Comm comm, const ufc_dofmap &dofmap, mesh::Topology &topology)
Create dof map on mesh from a ufc_dofmap.
Definition: utils.cpp:149
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:400
std::vector< std::vector< std::array< std::shared_ptr< const fem::FunctionSpace >, 2 > > > extract_function_spaces(const std::vector< std::vector< const fem::Form< T > * >> &a)
Extract test (0) and trial (1) function spaces pairs for each bilinear form for a rectangular array o...
Definition: utils.h:58
std::vector< std::string > get_coefficient_names(const ufc_form &ufc_form)
Get the name of each coefficient in a UFC form.
Definition: utils.cpp:179
std::vector< std::string > get_constant_names(const ufc_form &ufc_form)
Get the name of each constant in a UFC form.
Definition: utils.cpp:188
fem::CoordinateElement create_coordinate_map(const ufc_coordinate_mapping &ufc_cmap)
Create a CoordinateElement from ufc.
Definition: utils.cpp:198
IntegralType
Type of integral.
Definition: Form.h:27
la::SparsityPattern create_sparsity_pattern(const Form< T > &a)
Create a sparsity pattern for a given form. The pattern is not finalised, i.e. the caller is responsi...
Definition: utils.h:83
CellType
Cell type identifier.
Definition: cell_types.h:21