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 
56 template <typename T>
57 std::vector<
58  std::vector<std::array<std::shared_ptr<const function::FunctionSpace>, 2>>>
60  const Eigen::Ref<const Eigen::Array<const fem::Form<T>*, Eigen::Dynamic,
61  Eigen::Dynamic, Eigen::RowMajor>>& a)
62 {
63  std::vector<std::vector<
64  std::array<std::shared_ptr<const function::FunctionSpace>, 2>>>
65  spaces(a.rows(),
66  std::vector<
67  std::array<std::shared_ptr<const function::FunctionSpace>, 2>>(
68  a.cols()));
69  for (int i = 0; i < a.rows(); ++i)
70  for (int j = 0; j < a.cols(); ++j)
71  if (a(i, j))
72  spaces[i][j] = {a(i, j)->function_space(0), a(i, j)->function_space(1)};
73  return spaces;
74 }
75 
81 template <typename T>
83 {
84  if (a.rank() != 2)
85  {
86  throw std::runtime_error(
87  "Cannot create sparsity pattern. Form is not a bilinear form");
88  }
89 
90  // Get dof maps and mesh
91  std::array dofmaps{a.function_space(0)->dofmap().get(),
92  a.function_space(1)->dofmap().get()};
93  std::shared_ptr mesh = a.mesh();
94  assert(mesh);
95 
96  const std::set<IntegralType> types = a.integrals().types();
97  if (types.find(IntegralType::interior_facet) != types.end()
98  or types.find(IntegralType::exterior_facet) != types.end())
99  {
100  // FIXME: cleanup these calls? Some of the happen internally again.
101  const int tdim = mesh->topology().dim();
102  mesh->topology_mutable().create_entities(tdim - 1);
103  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
104  }
105 
106  return create_sparsity_pattern(mesh->topology(), dofmaps, types);
107 }
108 
114  const std::array<const DofMap*, 2>& dofmaps,
115  const std::set<IntegralType>& integrals);
116 
118 ElementDofLayout create_element_dof_layout(const ufc_dofmap& dofmap,
119  const mesh::CellType cell_type,
120  const std::vector<int>& parent_map
121  = {});
122 
127 DofMap create_dofmap(MPI_Comm comm, const ufc_dofmap& dofmap,
128  mesh::Topology& topology);
129 
131 template <typename T>
132 std::vector<
133  std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
134 get_coeffs_from_ufc_form(const ufc_form& ufc_form)
135 {
136  std::vector<
137  std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
138  coeffs;
139  const char** names = ufc_form.coefficient_name_map();
140  for (int i = 0; i < ufc_form.num_coefficients; ++i)
141  {
142  coeffs.emplace_back(ufc_form.original_coefficient_position(i), names[i],
143  nullptr);
144  }
145  return coeffs;
146 }
147 
149 template <typename T>
150 std::vector<
151  std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
152 get_constants_from_ufc_form(const ufc_form& ufc_form)
153 {
154  std::vector<
155  std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
156  constants;
157  const char** names = ufc_form.constant_name_map();
158  for (int i = 0; i < ufc_form.num_constants; ++i)
159  constants.emplace_back(names[i], nullptr);
160  return constants;
161 }
162 
166 template <typename T>
168  const ufc_form& ufc_form,
169  const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
170 {
171  assert(ufc_form.rank == (int)spaces.size());
172 
173  // Check argument function spaces
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 
188  // Get list of integral IDs, and load tabulate tensor into memory for each
189  bool needs_permutation_data = false;
190  std::map<IntegralType,
191  std::vector<std::pair<
192  int, std::function<void(T*, const T*, const T*, const double*,
193  const int*, const std::uint8_t*,
194  const std::uint32_t)>>>>
195  integral_data;
196 
197  std::vector<int> cell_integral_ids(ufc_form.num_cell_integrals);
198  ufc_form.get_cell_integral_ids(cell_integral_ids.data());
199  for (int id : cell_integral_ids)
200  {
201  ufc_integral* integral = ufc_form.create_cell_integral(id);
202  assert(integral);
203  if (integral->needs_permutation_data)
204  needs_permutation_data = true;
205  integral_data[IntegralType::cell].emplace_back(id,
206  integral->tabulate_tensor);
207  std::free(integral);
208  }
209 
210  // FIXME: Can this be handled better?
211  // FIXME: Handle forms with no space
212  if (ufc_form.num_exterior_facet_integrals > 0
213  or ufc_form.num_interior_facet_integrals > 0)
214  {
215  if (!spaces.empty())
216  {
217  auto mesh = spaces[0]->mesh();
218  const int tdim = mesh->topology().dim();
219  spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
220  }
221  }
222 
223  std::vector<int> exterior_facet_integral_ids(
224  ufc_form.num_exterior_facet_integrals);
225  ufc_form.get_exterior_facet_integral_ids(exterior_facet_integral_ids.data());
226  for (int id : exterior_facet_integral_ids)
227  {
228  ufc_integral* integral = ufc_form.create_exterior_facet_integral(id);
229  assert(integral);
230  if (integral->needs_permutation_data)
231  needs_permutation_data = true;
232  integral_data[IntegralType::exterior_facet].emplace_back(
233  id, integral->tabulate_tensor);
234  std::free(integral);
235  }
236 
237  std::vector<int> interior_facet_integral_ids(
238  ufc_form.num_interior_facet_integrals);
239  ufc_form.get_interior_facet_integral_ids(interior_facet_integral_ids.data());
240  for (int id : interior_facet_integral_ids)
241  {
242  ufc_integral* integral = ufc_form.create_interior_facet_integral(id);
243  assert(integral);
244  if (integral->needs_permutation_data)
245  needs_permutation_data = true;
246  integral_data[IntegralType::interior_facet].emplace_back(
247  id, integral->tabulate_tensor);
248  std::free(integral);
249  }
250 
251  // Not currently working
252  std::vector<int> vertex_integral_ids(ufc_form.num_vertex_integrals);
253  ufc_form.get_vertex_integral_ids(vertex_integral_ids.data());
254  if (vertex_integral_ids.size() > 0)
255  {
256  throw std::runtime_error(
257  "Vertex integrals not supported. Under development.");
258  }
259 
260  return fem::Form(
261  spaces, FormIntegrals<T>(integral_data, needs_permutation_data),
262  FormCoefficients<T>(fem::get_coeffs_from_ufc_form<T>(ufc_form)),
263  fem::get_constants_from_ufc_form<T>(ufc_form));
264 }
265 
272 template <typename T>
273 std::shared_ptr<Form<T>> create_form(
274  ufc_form* (*fptr)(),
275  const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
276 {
277  ufc_form* form = fptr();
278  auto L = std::make_shared<fem::Form<T>>(
279  dolfinx::fem::create_form<T>(*form, spaces));
280  std::free(form);
281  return L;
282 }
283 
288 create_coordinate_map(const ufc_coordinate_mapping& ufc_cmap);
289 
294 fem::CoordinateElement create_coordinate_map(ufc_coordinate_mapping* (*fptr)());
295 
305 std::shared_ptr<function::FunctionSpace>
306 create_functionspace(ufc_function_space* (*fptr)(const char*),
307  const std::string function_name,
308  std::shared_ptr<mesh::Mesh> mesh);
309 
310 // NOTE: This is subject to change
312 template <typename T>
313 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
315 {
316  // Get form coefficient offsets amd dofmaps
317  const fem::FormCoefficients<T>& coefficients = form.coefficients();
318  const std::vector<int>& offsets = coefficients.offsets();
319  std::vector<const fem::DofMap*> dofmaps(coefficients.size());
320  std::vector<Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>> v;
321  for (int i = 0; i < coefficients.size(); ++i)
322  {
323  dofmaps[i] = coefficients.get(i)->function_space()->dofmap().get();
324  v.emplace_back(coefficients.get(i)->x()->array());
325  }
326 
327  // Get mesh
328  std::shared_ptr<const mesh::Mesh> mesh = form.mesh();
329  assert(mesh);
330  const int tdim = mesh->topology().dim();
331  const std::int32_t num_cells
332  = mesh->topology().index_map(tdim)->size_local()
333  + mesh->topology().index_map(tdim)->num_ghosts();
334 
335  // Copy data into coefficient array
336  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
337  num_cells, offsets.back());
338  if (coefficients.size() > 0)
339  {
340  for (int cell = 0; cell < num_cells; ++cell)
341  {
342  for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
343  {
344  auto dofs = dofmaps[coeff]->cell_dofs(cell);
345  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& _v
346  = v[coeff];
347  for (Eigen::Index k = 0; k < dofs.size(); ++k)
348  c(cell, k + offsets[coeff]) = _v[dofs[k]];
349  }
350  }
351  }
352 
353  return c;
354 }
355 
356 // NOTE: This is subject to change
358 template <typename T>
359 Eigen::Array<T, Eigen::Dynamic, 1> pack_constants(const fem::Form<T>& form)
360 {
361  std::vector<T> constant_values;
362  for (auto& constant : form.constants())
363  {
364  const std::vector<T>& array = constant.second->value;
365  constant_values.insert(constant_values.end(), array.begin(), array.end());
366  }
367 
368  return Eigen::Map<const Eigen::Array<T, Eigen::Dynamic, 1>>(
369  constant_values.data(), constant_values.size(), 1);
370 }
371 
372 } // namespace fem
373 } // namespace dolfinx
dolfinx::fem::IntegralType
IntegralType
Type of integral.
Definition: FormIntegrals.h:27
dolfinx::fem::FormCoefficients::offsets
std::vector< int > offsets() const
Offset for each coefficient expansion array on a cell. Used to pack data for multiple coefficients in...
Definition: FormCoefficients.h:58
dolfinx::fem::FormCoefficients::size
int size() const
Get number of coefficients.
Definition: FormCoefficients.h:53
dolfinx::mesh::CellType
CellType
Cell type identifier.
Definition: cell_types.h:21
dolfinx::fem::Form::constants
std::vector< std::pair< std::string, std::shared_ptr< const function::Constant< T > > > > & constants()
Access constants.
Definition: Form.h:244
dolfinx::fem::Form::rank
int rank() const
Rank of the form (bilinear form = 2, linear form = 1, functional = 0, etc)
Definition: Form.h:129
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:59
dolfinx::fem::Form
Class for variational forms.
Definition: assembler.h:28
dolfinx::fem::create_coordinate_map
fem::CoordinateElement create_coordinate_map(const ufc_coordinate_mapping &ufc_cmap)
Create a CoordinateElement from ufc.
Definition: utils.cpp:210
dolfinx::fem::Form::coefficients
FormCoefficients< T > & coefficients()
Access coefficients.
Definition: Form.h:230
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:246
dolfinx::fem::Form::integrals
const FormIntegrals< T > & integrals() const
Access form integrals.
Definition: Form.h:236
dolfinx::fem::pack_constants
Eigen::Array< T, Eigen::Dynamic, 1 > pack_constants(const fem::Form< T > &form)
Pack form constants ready for assembly.
Definition: utils.h:359
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::function_space
std::shared_ptr< const function::FunctionSpace > function_space(int i) const
Return function space for given argument.
Definition: Form.h:204
dolfinx::fem::Form::mesh
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh from form.
Definition: Form.h:199
dolfinx::fem::pack_coefficients
Eigen::Array< T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > pack_coefficients(const fem::Form< T > &form)
Pack form coefficients ready for assembly.
Definition: utils.h:314
dolfinx::fem::get_constants_from_ufc_form
std::vector< std::pair< std::string, std::shared_ptr< const function::Constant< T > > > > get_constants_from_ufc_form(const ufc_form &ufc_form)
Extract coefficients from a UFC form.
Definition: utils.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:82
dolfinx::fem::FormCoefficients::get
std::shared_ptr< const function::Function< T > > get(int i) const
Get the Function coefficient i.
Definition: FormCoefficients.h:90
dolfinx::fem::get_coeffs_from_ufc_form
std::vector< std::tuple< int, std::string, std::shared_ptr< function::Function< T > > > > get_coeffs_from_ufc_form(const ufc_form &ufc_form)
Extract coefficients from a UFC form.
Definition: utils.h:134
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::FormCoefficients
Storage for the coefficients of a Form consisting of Function and the Element objects they are define...
Definition: FormCoefficients.h:34
dolfinx::fem::create_form
Form< T > create_form(const ufc_form &ufc_form, const std::vector< std::shared_ptr< const function::FunctionSpace >> &spaces)
Create a Form from UFC input.
Definition: utils.h:167
dolfinx::fem::CoordinateElement
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:24