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/function/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 function
32 {
33 template <typename T>
34 class Constant;
35 template <typename T>
36 class Function;
37 class FunctionSpace;
38 } // namespace function
39 
40 namespace mesh
41 {
42 class Mesh;
43 class Topology;
44 } // namespace mesh
45 
46 namespace fem
47 {
48 
49 template <typename T>
50 class Form;
51 
59 template <typename T>
60 std::vector<
61  std::vector<std::array<std::shared_ptr<const function::FunctionSpace>, 2>>>
63  const Eigen::Ref<const Eigen::Array<const fem::Form<T>*, Eigen::Dynamic,
64  Eigen::Dynamic, Eigen::RowMajor>>& a)
65 {
66  std::vector<std::vector<
67  std::array<std::shared_ptr<const function::FunctionSpace>, 2>>>
68  spaces(a.rows(),
69  std::vector<
70  std::array<std::shared_ptr<const function::FunctionSpace>, 2>>(
71  a.cols()));
72  for (int i = 0; i < a.rows(); ++i)
73  {
74  for (int j = 0; j < a.cols(); ++j)
75  {
76  if (a(i, j))
77  {
78  spaces[i][j]
79  = {a(i, j)->function_spaces()[0], a(i, j)->function_spaces()[1]};
80  }
81  }
82  }
83  return spaces;
84 }
85 
91 template <typename T>
93 {
94  if (a.rank() != 2)
95  {
96  throw std::runtime_error(
97  "Cannot create sparsity pattern. Form is not a bilinear form");
98  }
99 
100  // Get dof maps and mesh
101  std::array dofmaps{a.function_spaces().at(0)->dofmap().get(),
102  a.function_spaces().at(1)->dofmap().get()};
103  std::shared_ptr mesh = a.mesh();
104  assert(mesh);
105 
106  const std::set<IntegralType> types = a.integral_types();
107  if (types.find(IntegralType::interior_facet) != types.end()
108  or types.find(IntegralType::exterior_facet) != types.end())
109  {
110  // FIXME: cleanup these calls? Some of the happen internally again.
111  const int tdim = mesh->topology().dim();
112  mesh->topology_mutable().create_entities(tdim - 1);
113  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
114  }
115 
116  return create_sparsity_pattern(mesh->topology(), dofmaps, types);
117 }
118 
124  const std::array<const DofMap*, 2>& dofmaps,
125  const std::set<IntegralType>& integrals);
126 
128 ElementDofLayout create_element_dof_layout(const ufc_dofmap& dofmap,
129  const mesh::CellType cell_type,
130  const std::vector<int>& parent_map
131  = {});
132 
137 DofMap create_dofmap(MPI_Comm comm, const ufc_dofmap& dofmap,
138  mesh::Topology& topology);
139 
143 std::vector<std::string> get_coefficient_names(const ufc_form& ufc_form);
144 
148 std::vector<std::string> get_constant_names(const ufc_form& ufc_form);
149 
157 template <typename T>
159  const ufc_form& ufc_form,
160  const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces,
161  const std::vector<std::shared_ptr<const function::Function<T>>>&
162  coefficients,
163  const std::vector<std::shared_ptr<const function::Constant<T>>>& constants,
164  const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
165  const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
166 {
167  if (ufc_form.rank != (int)spaces.size())
168  throw std::runtime_error("Wrong number of argument spaces for Form.");
169  if (ufc_form.num_coefficients != (int)coefficients.size())
170  {
171  throw std::runtime_error(
172  "Mismatch between number of expected and provided Form coefficients.");
173  }
174  if (ufc_form.num_constants != (int)constants.size())
175  {
176  throw std::runtime_error(
177  "Mismatch between number of expected and provided Form constants.");
178  }
179 
180  // Check argument function spaces
181 #ifdef DEBUG
182  for (std::size_t i = 0; i < spaces.size(); ++i)
183  {
184  assert(spaces[i]->element());
185  std::unique_ptr<ufc_finite_element, decltype(free)*> ufc_element(
186  ufc_form.create_finite_element(i), free);
187  assert(ufc_element);
188  if (std::string(ufc_element->signature)
189  != spaces[i]->element()->signature())
190  {
191  throw std::runtime_error(
192  "Cannot create form. Wrong type of function space for argument.");
193  }
194  }
195 #endif
196 
197  // Get list of integral IDs, and load tabulate tensor into memory for each
198  using kern = std::function<void(PetscScalar*, const PetscScalar*,
199  const PetscScalar*, const double*, const int*,
200  const std::uint8_t*, const std::uint32_t)>;
201  std::map<IntegralType, std::pair<std::vector<std::pair<int, kern>>,
202  const mesh::MeshTags<int>*>>
203  integral_data;
204 
205  bool needs_permutation_data = false;
206 
207  // Attach cell kernels
208  std::vector<int> cell_integral_ids(ufc_form.num_cell_integrals);
209  ufc_form.get_cell_integral_ids(cell_integral_ids.data());
210  for (int id : cell_integral_ids)
211  {
212  ufc_integral* integral = ufc_form.create_cell_integral(id);
213  assert(integral);
214  if (integral->needs_permutation_data)
215  needs_permutation_data = true;
216  integral_data[IntegralType::cell].first.emplace_back(
217  id, integral->tabulate_tensor);
218  std::free(integral);
219  }
220 
221  // Attach cell subdomain data
222  if (auto it = subdomains.find(IntegralType::cell);
223  it != subdomains.end() and !cell_integral_ids.empty())
224  {
225  integral_data[IntegralType::cell].second = it->second;
226  }
227 
228  // FIXME: Can facets be handled better?
229 
230  // Create facets, if required
231  if (ufc_form.num_exterior_facet_integrals > 0
232  or ufc_form.num_interior_facet_integrals > 0)
233  {
234  if (!spaces.empty())
235  {
236  auto mesh = spaces[0]->mesh();
237  const int tdim = mesh->topology().dim();
238  spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
239  }
240  }
241 
242  // Attach exterior facet kernels
243  std::vector<int> exterior_facet_integral_ids(
244  ufc_form.num_exterior_facet_integrals);
245  ufc_form.get_exterior_facet_integral_ids(exterior_facet_integral_ids.data());
246  for (int id : exterior_facet_integral_ids)
247  {
248  ufc_integral* integral = ufc_form.create_exterior_facet_integral(id);
249  assert(integral);
250  if (integral->needs_permutation_data)
251  needs_permutation_data = true;
252  integral_data[IntegralType::exterior_facet].first.emplace_back(
253  id, integral->tabulate_tensor);
254  std::free(integral);
255  }
256 
257  // Attach exterior facet subdomain data
258  if (auto it = subdomains.find(IntegralType::exterior_facet);
259  it != subdomains.end() and !exterior_facet_integral_ids.empty())
260  {
261  integral_data[IntegralType::exterior_facet].second = it->second;
262  }
263 
264  // Attach interior facet kernels
265  std::vector<int> interior_facet_integral_ids(
266  ufc_form.num_interior_facet_integrals);
267  ufc_form.get_interior_facet_integral_ids(interior_facet_integral_ids.data());
268  for (int id : interior_facet_integral_ids)
269  {
270  ufc_integral* integral = ufc_form.create_interior_facet_integral(id);
271  assert(integral);
272  if (integral->needs_permutation_data)
273  needs_permutation_data = true;
274  integral_data[IntegralType::interior_facet].first.emplace_back(
275  id, integral->tabulate_tensor);
276  std::free(integral);
277  }
278 
279  // Attach interior facet subdomain data
280  if (auto it = subdomains.find(IntegralType::interior_facet);
281  it != subdomains.end() and !interior_facet_integral_ids.empty())
282  {
283  integral_data[IntegralType::interior_facet].second = it->second;
284  }
285 
286  // Vertex integrals: not currently working
287  std::vector<int> vertex_integral_ids(ufc_form.num_vertex_integrals);
288  ufc_form.get_vertex_integral_ids(vertex_integral_ids.data());
289  if (!vertex_integral_ids.empty())
290  {
291  throw std::runtime_error(
292  "Vertex integrals not supported. Under development.");
293  }
294 
295  return fem::Form(spaces, integral_data, coefficients, constants,
296  needs_permutation_data, mesh);
297 }
298 
308 template <typename T>
310  const ufc_form& ufc_form,
311  const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces,
312  const std::map<std::string, std::shared_ptr<const function::Function<T>>>&
313  coefficients,
314  const std::map<std::string, std::shared_ptr<const function::Constant<T>>>&
315  constants,
316  const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
317  const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
318 {
319  // Place coefficients in appropriate order
320  std::vector<std::shared_ptr<const function::Function<T>>> coeff_map;
321  for (const std::string& name : get_coefficient_names(ufc_form))
322  {
323  if (auto it = coefficients.find(name); it != coefficients.end())
324  coeff_map.push_back(it->second);
325  else
326  {
327  throw std::runtime_error("Form coefficient \"" + name
328  + "\" not provided.");
329  }
330  }
331 
332  // Place constants in appropriate order
333  std::vector<std::shared_ptr<const function::Constant<T>>> const_map;
334  for (const std::string& name : get_constant_names(ufc_form))
335  {
336  if (auto it = constants.find(name); it != constants.end())
337  const_map.push_back(it->second);
338  else
339  {
340  throw std::runtime_error("Form constant \"" + name + "\" not provided.");
341  }
342  }
343 
344  return create_form(ufc_form, spaces, coeff_map, const_map, subdomains, mesh);
345 }
346 
358 template <typename T>
359 std::shared_ptr<Form<T>> create_form(
360  ufc_form* (*fptr)(),
361  const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces,
362  const std::map<std::string, std::shared_ptr<const function::Function<T>>>&
363  coefficients,
364  const std::map<std::string, std::shared_ptr<const function::Constant<T>>>&
365  constants,
366  const std::map<IntegralType, const mesh::MeshTags<int>*>& subdomains,
367  const std::shared_ptr<const mesh::Mesh>& mesh = nullptr)
368 {
369  ufc_form* form = fptr();
370  auto L = std::make_shared<fem::Form<T>>(dolfinx::fem::create_form<T>(
371  *form, spaces, coefficients, constants, subdomains, mesh));
372  std::free(form);
373  return L;
374 }
375 
380 create_coordinate_map(const ufc_coordinate_mapping& ufc_cmap);
381 
386 fem::CoordinateElement create_coordinate_map(ufc_coordinate_mapping* (*fptr)());
387 
397 std::shared_ptr<function::FunctionSpace>
398 create_functionspace(ufc_function_space* (*fptr)(const char*),
399  const std::string function_name,
400  std::shared_ptr<mesh::Mesh> mesh);
401 
402 // NOTE: This is subject to change
404 template <typename U>
405 Eigen::Array<typename U::scalar_type, Eigen::Dynamic, Eigen::Dynamic,
406  Eigen::RowMajor>
407 pack_coefficients(const U& u)
408 {
409  using T = typename U::scalar_type;
410 
411  // Get form coefficient offsets and dofmaps
412  const std::vector<std::shared_ptr<const function::Function<T>>> coefficients
413  = u.coefficients();
414  const std::vector<int> offsets = u.coefficient_offsets();
415  std::vector<const fem::DofMap*> dofmaps(coefficients.size());
416  std::vector<Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>> v;
417  for (std::size_t i = 0; i < coefficients.size(); ++i)
418  {
419  dofmaps[i] = coefficients[i]->function_space()->dofmap().get();
420  v.emplace_back(coefficients[i]->x()->array());
421  }
422 
423  // Get mesh
424  std::shared_ptr<const mesh::Mesh> mesh = u.mesh();
425  assert(mesh);
426  const int tdim = mesh->topology().dim();
427  const std::int32_t num_cells
428  = mesh->topology().index_map(tdim)->size_local()
429  + mesh->topology().index_map(tdim)->num_ghosts();
430 
431  // Copy data into coefficient array
432  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
433  num_cells, offsets.back());
434  if (coefficients.size() > 0)
435  {
436  for (int cell = 0; cell < num_cells; ++cell)
437  {
438  for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
439  {
440  auto dofs = dofmaps[coeff]->cell_dofs(cell);
441  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& _v
442  = v[coeff];
443  for (Eigen::Index k = 0; k < dofs.size(); ++k)
444  c(cell, k + offsets[coeff]) = _v[dofs[k]];
445  }
446  }
447  }
448 
449  return c;
450 }
451 
452 // NOTE: This is subject to change
454 template <typename U>
455 Eigen::Array<typename U::scalar_type, Eigen::Dynamic, 1>
456 pack_constants(const U& u)
457 {
458  using T = typename U::scalar_type;
459 
460  const auto& constants = u.constants();
461 
462  // Calculate size of array needed to store packed constants.
463  Eigen::Index size
464  = std::accumulate(constants.begin(), constants.end(), 0,
465  [](Eigen::Index sum, const auto& constant) {
466  return sum + constant->value.size();
467  });
468 
469  // Pack constants.
470  Eigen::Array<T, Eigen::Dynamic, 1> constant_values(size);
471  Eigen::Index offset = 0;
472  for (const auto& constant : constants)
473  {
474  const auto& value = constant->value;
475  const Eigen::Index value_size = value.size();
476 
477  for (Eigen::Index i = 0; i < value_size; ++i)
478  constant_values[offset + i] = value[i];
479 
480  offset += value_size;
481  }
482 
483  return constant_values;
484 }
485 
486 } // namespace fem
487 } // namespace dolfinx
dolfinx::fem::pack_coefficients
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:407
dolfinx::fem::IntegralType
IntegralType
Type of integral.
Definition: Form.h:33
dolfinx::mesh::CellType
CellType
Cell type identifier.
Definition: cell_types.h:21
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::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::fem::get_constant_names
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:218
dolfinx::fem::get_coefficient_names
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:209
dolfinx::la::SparsityPattern
This class provides a sparsity pattern data structure that can be used to initialize sparse matrices.
Definition: SparsityPattern.h:37
dolfinx::fem::extract_function_spaces
std::vector< std::vector< std::array< std::shared_ptr< const function::FunctionSpace >, 2 > > > extract_function_spaces(const Eigen::Ref< const Eigen::Array< const fem::Form< T > *, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> &a)
Extract test (0) and trial (1) function spaces pairs for each bilinear form for a rectangular array o...
Definition: utils.h:62
dolfinx::fem::Form
Class for variational forms.
Definition: Form.h:66
dolfinx::fem::create_coordinate_map
fem::CoordinateElement create_coordinate_map(const ufc_coordinate_mapping &ufc_cmap)
Create a CoordinateElement from ufc.
Definition: utils.cpp:228
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::function::Constant
A constant value which can be attached to a Form. Constants may be scalar (rank 0),...
Definition: Constant.h:19
dolfinx::fem::create_functionspace
std::shared_ptr< function::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:264
dolfinx::fem::create_form
Form< T > create_form(const ufc_form &ufc_form, const std::vector< std::shared_ptr< const function::FunctionSpace >> &spaces, const std::vector< std::shared_ptr< const function::Function< T >>> &coefficients, const std::vector< std::shared_ptr< const function::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:158
dolfinx::fem::create_element_dof_layout
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:100
dolfinx::fem::create_dofmap
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:179
dolfinx::fem::Form::mesh
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh for the form.
Definition: Form.h:152
dolfinx::fem::create_sparsity_pattern
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:92
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::pack_constants
Eigen::Array< typename U::scalar_type, Eigen::Dynamic, 1 > pack_constants(const U &u)
Pack constants of u of generic type U ready for assembly.
Definition: utils.h:456
dolfinx::fem::CoordinateElement
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:24