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 <utility>
21 #include <vector>
22 
23 struct ufc_dofmap;
24 struct ufc_form;
25 struct ufc_coordinate_mapping;
26 struct ufc_function_space;
27 
28 namespace dolfinx
29 {
30 namespace common
31 {
32 class IndexMap;
33 }
34 
35 namespace function
36 {
37 template <typename T>
38 class Constant;
39 template <typename T>
40 class Function;
41 class FunctionSpace;
42 } // namespace function
43 
44 namespace mesh
45 {
46 class Mesh;
47 class Topology;
48 } // namespace mesh
49 
50 namespace fem
51 {
52 
60 template <typename T>
61 Eigen::Array<std::array<std::shared_ptr<const function::FunctionSpace>, 2>,
62  Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
64  const Eigen::Ref<const Eigen::Array<const fem::Form<T>*, Eigen::Dynamic,
65  Eigen::Dynamic, Eigen::RowMajor>>& a)
66 {
67  Eigen::Array<std::array<std::shared_ptr<const function::FunctionSpace>, 2>,
68  Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
69  spaces(a.rows(), a.cols());
70  for (int i = 0; i < a.rows(); ++i)
71  for (int j = 0; j < a.cols(); ++j)
72  if (a(i, j))
73  spaces(i, j) = {a(i, j)->function_space(0), a(i, j)->function_space(1)};
74  return spaces;
75 }
76 
86 std::array<std::vector<std::shared_ptr<const function::FunctionSpace>>, 2>
88  const Eigen ::Ref<const Eigen::Array<
89  std::array<std::shared_ptr<const function::FunctionSpace>, 2>,
90  Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>>& V);
91 
97 template <typename T>
99 {
100  if (a.rank() != 2)
101  {
102  throw std::runtime_error(
103  "Cannot create sparsity pattern. Form is not a bilinear form");
104  }
105 
106  // Get dof maps and mesh
107  std::array dofmaps{a.function_space(0)->dofmap().get(),
108  a.function_space(1)->dofmap().get()};
109  std::shared_ptr mesh = a.mesh();
110  assert(mesh);
111 
112  const std::set<IntegralType> types = a.integrals().types();
113  if (types.find(IntegralType::interior_facet) != types.end()
114  or types.find(IntegralType::exterior_facet) != types.end())
115  {
116  // FIXME: cleanup these calls? Some of the happen internally again.
117  const int tdim = mesh->topology().dim();
118  mesh->topology_mutable().create_entities(tdim - 1);
119  mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
120  }
121 
122  return create_sparsity_pattern(mesh->topology(), dofmaps, types);
123 }
124 
130  const std::array<const DofMap*, 2>& dofmaps,
131  const std::set<IntegralType>& integrals);
132 
134 ElementDofLayout create_element_dof_layout(const ufc_dofmap& dofmap,
135  const mesh::CellType cell_type,
136  const std::vector<int>& parent_map
137  = {});
138 
143 DofMap create_dofmap(MPI_Comm comm, const ufc_dofmap& dofmap,
144  mesh::Topology& topology);
145 
147 template <typename T>
148 std::vector<
149  std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
150 get_coeffs_from_ufc_form(const ufc_form& ufc_form)
151 {
152  std::vector<
153  std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
154  coeffs;
155  const char** names = ufc_form.coefficient_name_map();
156  for (int i = 0; i < ufc_form.num_coefficients; ++i)
157  {
158  coeffs.emplace_back(ufc_form.original_coefficient_position(i), names[i],
159  nullptr);
160  }
161  return coeffs;
162 }
163 
165 template <typename T>
166 std::vector<
167  std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
168 get_constants_from_ufc_form(const ufc_form& ufc_form)
169 {
170  std::vector<
171  std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
172  constants;
173  const char** names = ufc_form.constant_name_map();
174  for (int i = 0; i < ufc_form.num_constants; ++i)
175  constants.emplace_back(names[i], nullptr);
176  return constants;
177 }
178 
182 template <typename T>
184  const ufc_form& ufc_form,
185  const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
186 {
187  assert(ufc_form.rank == (int)spaces.size());
188 
189  // Check argument function spaces
190  for (std::size_t i = 0; i < spaces.size(); ++i)
191  {
192  assert(spaces[i]->element());
193  std::unique_ptr<ufc_finite_element, decltype(free)*> ufc_element(
194  ufc_form.create_finite_element(i), free);
195  assert(ufc_element);
196  if (std::string(ufc_element->signature)
197  != spaces[i]->element()->signature())
198  {
199  throw std::runtime_error(
200  "Cannot create form. Wrong type of function space for argument.");
201  }
202  }
203 
204  // Get list of integral IDs, and load tabulate tensor into memory for each
205  FormIntegrals<T> integrals;
206 
207  std::vector<int> cell_integral_ids(ufc_form.num_cell_integrals);
208  ufc_form.get_cell_integral_ids(cell_integral_ids.data());
209  for (int id : cell_integral_ids)
210  {
211  ufc_integral* cell_integral = ufc_form.create_cell_integral(id);
212  assert(cell_integral);
213  integrals.set_tabulate_tensor(IntegralType::cell, id,
214  cell_integral->tabulate_tensor);
215  std::free(cell_integral);
216  }
217 
218  // FIXME: Can this be handled better?
219  // FIXME: Handle forms with no space
220  if (ufc_form.num_exterior_facet_integrals > 0
221  or ufc_form.num_interior_facet_integrals > 0)
222  {
223  if (!spaces.empty())
224  {
225  auto mesh = spaces[0]->mesh();
226  const int tdim = mesh->topology().dim();
227  spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
228  }
229  }
230 
231  std::vector<int> exterior_facet_integral_ids(
232  ufc_form.num_exterior_facet_integrals);
233  ufc_form.get_exterior_facet_integral_ids(exterior_facet_integral_ids.data());
234  for (int id : exterior_facet_integral_ids)
235  {
236  ufc_integral* exterior_facet_integral
237  = ufc_form.create_exterior_facet_integral(id);
238  assert(exterior_facet_integral);
239  integrals.set_tabulate_tensor(IntegralType::exterior_facet, id,
240  exterior_facet_integral->tabulate_tensor);
241  std::free(exterior_facet_integral);
242  }
243 
244  std::vector<int> interior_facet_integral_ids(
245  ufc_form.num_interior_facet_integrals);
246  ufc_form.get_interior_facet_integral_ids(interior_facet_integral_ids.data());
247  for (int id : interior_facet_integral_ids)
248  {
249  ufc_integral* interior_facet_integral
250  = ufc_form.create_interior_facet_integral(id);
251  assert(interior_facet_integral);
252  integrals.set_tabulate_tensor(IntegralType::interior_facet, id,
253  interior_facet_integral->tabulate_tensor);
254  std::free(interior_facet_integral);
255  }
256 
257  // Not currently working
258  std::vector<int> vertex_integral_ids(ufc_form.num_vertex_integrals);
259  ufc_form.get_vertex_integral_ids(vertex_integral_ids.data());
260  if (vertex_integral_ids.size() > 0)
261  {
262  throw std::runtime_error(
263  "Vertex integrals not supported. Under development.");
264  }
265 
266  return fem::Form(
267  spaces, integrals,
268  FormCoefficients<T>(fem::get_coeffs_from_ufc_form<T>(ufc_form)),
269  fem::get_constants_from_ufc_form<T>(ufc_form));
270 }
271 
278 template <typename T>
279 std::shared_ptr<Form<T>> create_form(
280  ufc_form* (*fptr)(),
281  const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
282 {
283  ufc_form* form = fptr();
284  auto L = std::make_shared<fem::Form<T>>(
285  dolfinx::fem::create_form<T>(*form, spaces));
286  std::free(form);
287  return L;
288 }
289 
294 create_coordinate_map(const ufc_coordinate_mapping& ufc_cmap);
295 
300 fem::CoordinateElement create_coordinate_map(ufc_coordinate_mapping* (*fptr)());
301 
311 std::shared_ptr<function::FunctionSpace>
312 create_functionspace(ufc_function_space* (*fptr)(const char*),
313  const std::string function_name,
314  std::shared_ptr<mesh::Mesh> mesh);
315 
316 // NOTE: This is subject to change
318 template <typename T>
319 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
321 {
322  // Get form coefficient offsets amd dofmaps
323  const fem::FormCoefficients<T>& coefficients = form.coefficients();
324  const std::vector<int>& offsets = coefficients.offsets();
325  std::vector<const fem::DofMap*> dofmaps(coefficients.size());
326  std::vector<Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>> v;
327  for (int i = 0; i < coefficients.size(); ++i)
328  {
329  dofmaps[i] = coefficients.get(i)->function_space()->dofmap().get();
330  v.emplace_back(coefficients.get(i)->x()->array());
331  }
332 
333  // Get mesh
334  std::shared_ptr<const mesh::Mesh> mesh = form.mesh();
335  assert(mesh);
336  const int tdim = mesh->topology().dim();
337  const std::int32_t num_cells
338  = mesh->topology().index_map(tdim)->size_local()
339  + mesh->topology().index_map(tdim)->num_ghosts();
340 
341  // Copy data into coefficient array
342  Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
343  num_cells, offsets.back());
344  if (coefficients.size() > 0)
345  {
346  for (int cell = 0; cell < num_cells; ++cell)
347  {
348  for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
349  {
350  auto dofs = dofmaps[coeff]->cell_dofs(cell);
351  const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& _v
352  = v[coeff];
353  for (Eigen::Index k = 0; k < dofs.size(); ++k)
354  c(cell, k + offsets[coeff]) = _v[dofs[k]];
355  }
356  }
357  }
358 
359  return c;
360 }
361 
362 // NOTE: This is subject to change
364 template <typename T>
365 Eigen::Array<T, Eigen::Dynamic, 1> pack_constants(const fem::Form<T>& form)
366 {
367  std::vector<T> constant_values;
368  for (auto& constant : form.constants())
369  {
370  const std::vector<T>& array = constant.second->value;
371  constant_values.insert(constant_values.end(), array.begin(), array.end());
372  }
373 
374  return Eigen::Map<const Eigen::Array<T, Eigen::Dynamic, 1>>(
375  constant_values.data(), constant_values.size(), 1);
376 }
377 
378 } // namespace fem
379 } // namespace dolfinx
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:22
dolfinx::fem::Form::rank
int rank() const
Rank of form (bilinear form = 2, linear form = 1, functional = 0, etc)
Definition: Form.h:123
dolfinx::la::SparsityPattern
This class provides a sparsity pattern data structure that can be used to initialize sparse matrices.
Definition: SparsityPattern.h:36
dolfinx::fem::Form
Class for variational forms.
Definition: assembler.h:29
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::fem::create_coordinate_map
fem::CoordinateElement create_coordinate_map(const ufc_coordinate_mapping &ufc_cmap)
Create a CoordinateElement from ufc.
Definition: utils.cpp:291
dolfinx::fem::Form::coefficients
FormCoefficients< T > & coefficients()
Access coefficients.
Definition: Form.h:261
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:327
dolfinx::fem::Form::constants
const std::vector< std::pair< std::string, std::shared_ptr< const function::Constant< T > > > > & constants() const
Access constants.
Definition: Form.h:275
dolfinx::fem::Form::integrals
const FormIntegrals< T > & integrals() const
Access form integrals.
Definition: Form.h:267
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:365
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:176
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:251
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:235
dolfinx::fem::Form::mesh
std::shared_ptr< const mesh::Mesh > mesh() const
Extract common mesh from form.
Definition: Form.h:230
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:320
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:168
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:98
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:150
dolfinx::mesh::Topology
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition: Topology.h:57
dolfinx::fem::FormIntegrals
Integrals of a Form, including those defined over cells, interior and exterior facets,...
Definition: FormIntegrals.h:38
dolfinx::fem::FormCoefficients
Storage for the coefficients of a Form consisting of Function and the Element objects they are define...
Definition: FormCoefficients.h:33
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:183
dolfinx::fem::extract_function_spaces
Eigen::Array< std::array< std::shared_ptr< const function::FunctionSpace >, 2 >, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > 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:63
dolfinx::fem::common_function_spaces
std::array< std::vector< std::shared_ptr< const function::FunctionSpace > >, 2 > common_function_spaces(const Eigen ::Ref< const Eigen::Array< std::array< std::shared_ptr< const function::FunctionSpace >, 2 >, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor >> &V)
Extract FunctionSpaces for (0) rows blocks and (1) columns blocks from a rectangular array of bilinea...
Definition: utils.cpp:91
dolfinx::fem::CoordinateElement
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:23