9 #include "CoordinateElement.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>
57 std::vector<std::array<std::shared_ptr<const fem::FunctionSpace>, 2>>>
61 std::vector<std::array<std::shared_ptr<const fem::FunctionSpace>, 2>>>
64 std::vector<std::array<std::shared_ptr<const fem::FunctionSpace>, 2>>(
66 for (std::size_t i = 0; i < a.size(); ++i)
68 for (std::size_t j = 0; j < a[i].size(); ++j)
71 spaces[i][j] = {form->function_spaces()[0], form->function_spaces()[1]};
87 throw std::runtime_error(
88 "Cannot create sparsity pattern. Form is not a bilinear form");
92 std::array<const std::reference_wrapper<const fem::DofMap>, 2> dofmaps{
95 std::shared_ptr mesh = a.
mesh();
99 if (types.find(IntegralType::interior_facet) != types.end()
100 or types.find(IntegralType::exterior_facet) != types.end())
103 const int tdim = mesh->topology().dim();
104 mesh->topology_mutable().create_entities(tdim - 1);
105 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
116 const std::array<
const std::reference_wrapper<const fem::DofMap>, 2>&
118 const std::set<IntegralType>& integrals);
123 const std::vector<int>& parent_map
130 DofMap
create_dofmap(MPI_Comm comm,
const ufc_dofmap& dofmap,
131 mesh::Topology& topology);
150 template <
typename T>
152 const ufc_form& ufc_form,
153 const std::vector<std::shared_ptr<const fem::FunctionSpace>>& spaces,
157 const std::shared_ptr<const mesh::Mesh>& mesh =
nullptr)
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())
163 throw std::runtime_error(
164 "Mismatch between number of expected and provided Form coefficients.");
166 if (ufc_form.num_constants != (
int)constants.size())
168 throw std::runtime_error(
169 "Mismatch between number of expected and provided Form constants.");
174 for (std::size_t i = 0; i < spaces.size(); ++i)
176 assert(spaces[i]->element());
177 std::unique_ptr<ufc_finite_element, decltype(free)*> ufc_element(
178 ufc_form.create_finite_element(i), free);
180 if (std::string(ufc_element->signature)
181 != spaces[i]->element()->signature())
183 throw std::runtime_error(
184 "Cannot create form. Wrong type of function space for argument.");
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>>,
198 bool needs_permutation_data =
false;
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)
205 ufc_integral* integral = ufc_form.create_cell_integral(
id);
207 if (integral->needs_transformation_data)
208 needs_permutation_data =
true;
209 integral_data[IntegralType::cell].first.emplace_back(
210 id, integral->tabulate_tensor);
215 if (
auto it = subdomains.find(IntegralType::cell);
216 it != subdomains.end() and !cell_integral_ids.empty())
218 integral_data[IntegralType::cell].second = it->second;
224 if (ufc_form.num_exterior_facet_integrals > 0
225 or ufc_form.num_interior_facet_integrals > 0)
229 auto mesh = spaces[0]->mesh();
230 const int tdim = mesh->topology().dim();
231 spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
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)
241 ufc_integral* integral = ufc_form.create_exterior_facet_integral(
id);
243 if (integral->needs_transformation_data)
244 needs_permutation_data =
true;
245 integral_data[IntegralType::exterior_facet].first.emplace_back(
246 id, integral->tabulate_tensor);
251 if (
auto it = subdomains.find(IntegralType::exterior_facet);
252 it != subdomains.end() and !exterior_facet_integral_ids.empty())
254 integral_data[IntegralType::exterior_facet].second = it->second;
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)
263 ufc_integral* integral = ufc_form.create_interior_facet_integral(
id);
265 if (integral->needs_transformation_data)
266 needs_permutation_data =
true;
267 integral_data[IntegralType::interior_facet].first.emplace_back(
268 id, integral->tabulate_tensor);
273 if (
auto it = subdomains.find(IntegralType::interior_facet);
274 it != subdomains.end() and !interior_facet_integral_ids.empty())
276 integral_data[IntegralType::interior_facet].second = it->second;
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())
284 throw std::runtime_error(
285 "Vertex integrals not supported. Under development.");
288 return fem::Form(spaces, integral_data, coefficients, constants,
289 needs_permutation_data, mesh);
301 template <
typename T>
303 const ufc_form& ufc_form,
304 const std::vector<std::shared_ptr<const fem::FunctionSpace>>& spaces,
310 const std::shared_ptr<const mesh::Mesh>& mesh =
nullptr)
313 std::vector<std::shared_ptr<const fem::Function<T>>> coeff_map;
316 if (
auto it = coefficients.find(name); it != coefficients.end())
317 coeff_map.push_back(it->second);
320 throw std::runtime_error(
"Form coefficient \"" + name
321 +
"\" not provided.");
326 std::vector<std::shared_ptr<const fem::Constant<T>>> const_map;
329 if (
auto it = constants.find(name); it != constants.end())
330 const_map.push_back(it->second);
333 throw std::runtime_error(
"Form constant \"" + name +
"\" not provided.");
337 return create_form(ufc_form, spaces, coeff_map, const_map, subdomains, mesh);
351 template <
typename T>
354 const std::vector<std::shared_ptr<const fem::FunctionSpace>>& spaces,
360 const std::shared_ptr<const mesh::Mesh>& mesh =
nullptr)
362 ufc_form* form = fptr();
363 auto L = std::make_shared<fem::Form<T>>(fem::create_form<T>(
364 *form, spaces, coefficients, constants, subdomains, mesh));
390 std::shared_ptr<fem::FunctionSpace>
392 const std::string function_name,
393 std::shared_ptr<mesh::Mesh> mesh);
397 template <
typename U>
400 using T =
typename U::scalar_type;
403 const std::vector<std::shared_ptr<const fem::Function<T>>> coefficients
405 const std::vector<int> offsets = u.coefficient_offsets();
406 std::vector<const fem::DofMap*> dofmaps(coefficients.size());
407 std::vector<int> bs(coefficients.size());
408 std::vector<std::reference_wrapper<const std::vector<T>>> v;
409 v.reserve(coefficients.size());
410 for (std::size_t i = 0; i < coefficients.size(); ++i)
412 dofmaps[i] = coefficients[i]->function_space()->dofmap().get();
413 bs[i] = dofmaps[i]->bs();
414 v.push_back(coefficients[i]->x()->array());
418 std::shared_ptr<const mesh::Mesh> mesh = u.mesh();
420 const int tdim = mesh->topology().dim();
421 const std::int32_t num_cells
422 = mesh->topology().index_map(tdim)->size_local()
423 + mesh->topology().index_map(tdim)->num_ghosts();
427 if (!coefficients.empty())
429 for (
int cell = 0; cell < num_cells; ++cell)
431 for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
433 tcb::span<const std::int32_t> dofs = dofmaps[coeff]->cell_dofs(cell);
434 const std::vector<T>& _v = v[coeff];
435 for (std::size_t i = 0; i < dofs.size(); ++i)
437 for (
int k = 0; k < bs[coeff]; ++k)
439 c(cell, bs[coeff] * i + k + offsets[coeff])
440 = _v[bs[coeff] * dofs[i] + k];
452 template <
typename U>
455 using T =
typename U::scalar_type;
456 const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants
461 = std::accumulate(constants.begin(), constants.end(), 0,
462 [](std::int32_t sum,
const auto& constant) {
463 return sum + constant->value.size();
467 std::vector<T> constant_values(size);
468 std::int32_t offset = 0;
469 for (
const auto& constant : constants)
471 const std::vector<T>& value = constant->value;
472 for (std::size_t i = 0; i < value.size(); ++i)
473 constant_values[offset + i] = value[i];
474 offset += value.size();
477 return constant_values;
This class provides a dynamic 2-dimensional row-wise array data structure.
Definition: array2d.h:21
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:31
This class represents a function in a finite element function space , given by.
Definition: Function.h:46
This class provides a sparsity pattern data structure that can be used to initialize sparse matrices.
Definition: SparsityPattern.h:36
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:453
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:244
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
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
array2d< typename U::scalar_type > pack_coefficients(const U &u)
Pack coefficients of u of generic type U ready for assembly.
Definition: utils.h:398
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:23