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_permutation_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_permutation_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_permutation_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>>(dolfinx::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>
398 Eigen::Array<
typename U::scalar_type, Eigen::Dynamic, Eigen::Dynamic,
402 using T =
typename U::scalar_type;
405 const std::vector<std::shared_ptr<const fem::Function<T>>> 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)
413 dofmaps[i] = coefficients[i]->function_space()->dofmap().get();
414 bs[i] = dofmaps[i]->bs();
415 v.push_back(coefficients[i]->x()->array());
419 std::shared_ptr<const mesh::Mesh> mesh = u.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();
427 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
428 num_cells, offsets.back());
429 if (coefficients.size() > 0)
431 for (
int cell = 0; cell < num_cells; ++cell)
433 for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
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)
439 for (
int k = 0; k < bs[coeff]; ++k)
441 c(cell, bs[coeff] * i + k + offsets[coeff])
442 = _v[bs[coeff] * dofs[i] + k];
454 template <
typename U>
457 using T =
typename U::scalar_type;
458 const std::vector<std::shared_ptr<const fem::Constant<T>>>& constants
463 = std::accumulate(constants.begin(), constants.end(), 0,
464 [](std::int32_t sum,
const auto& constant) {
465 return sum + constant->value.size();
469 std::vector<T> constant_values(size);
470 std::int32_t offset = 0;
471 for (
const auto& constant : constants)
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();
479 return constant_values;
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
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
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