9 #include "CoordinateElement.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>
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)
66 std::vector<std::vector<
67 std::array<std::shared_ptr<const function::FunctionSpace>, 2>>>
70 std::array<std::shared_ptr<const function::FunctionSpace>, 2>>(
72 for (
int i = 0; i < a.rows(); ++i)
74 for (
int j = 0; j < a.cols(); ++j)
79 = {a(i, j)->function_spaces()[0], a(i, j)->function_spaces()[1]};
96 throw std::runtime_error(
97 "Cannot create sparsity pattern. Form is not a bilinear form");
103 std::shared_ptr mesh = a.
mesh();
107 if (types.find(IntegralType::interior_facet) != types.end()
108 or types.find(IntegralType::exterior_facet) != types.end())
111 const int tdim = mesh->topology().dim();
112 mesh->topology_mutable().create_entities(tdim - 1);
113 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
124 const std::array<const DofMap*, 2>& dofmaps,
125 const std::set<IntegralType>& integrals);
130 const std::vector<int>& parent_map
137 DofMap
create_dofmap(MPI_Comm comm,
const ufc_dofmap& dofmap,
138 mesh::Topology& topology);
157 template <
typename T>
159 const ufc_form& ufc_form,
160 const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces,
165 const std::shared_ptr<const mesh::Mesh>& mesh =
nullptr)
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())
171 throw std::runtime_error(
172 "Mismatch between number of expected and provided Form coefficients.");
174 if (ufc_form.num_constants != (
int)constants.size())
176 throw std::runtime_error(
177 "Mismatch between number of expected and provided Form constants.");
182 for (std::size_t i = 0; i < spaces.size(); ++i)
184 assert(spaces[i]->element());
185 std::unique_ptr<ufc_finite_element, decltype(free)*> ufc_element(
186 ufc_form.create_finite_element(i), free);
188 if (std::string(ufc_element->signature)
189 != spaces[i]->element()->signature())
191 throw std::runtime_error(
192 "Cannot create form. Wrong type of function space for argument.");
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>>,
205 bool needs_permutation_data =
false;
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)
212 ufc_integral* integral = ufc_form.create_cell_integral(
id);
214 if (integral->needs_permutation_data)
215 needs_permutation_data =
true;
216 integral_data[IntegralType::cell].first.emplace_back(
217 id, integral->tabulate_tensor);
222 if (
auto it = subdomains.find(IntegralType::cell);
223 it != subdomains.end() and !cell_integral_ids.empty())
225 integral_data[IntegralType::cell].second = it->second;
231 if (ufc_form.num_exterior_facet_integrals > 0
232 or ufc_form.num_interior_facet_integrals > 0)
236 auto mesh = spaces[0]->mesh();
237 const int tdim = mesh->topology().dim();
238 spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
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)
248 ufc_integral* integral = ufc_form.create_exterior_facet_integral(
id);
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);
258 if (
auto it = subdomains.find(IntegralType::exterior_facet);
259 it != subdomains.end() and !exterior_facet_integral_ids.empty())
261 integral_data[IntegralType::exterior_facet].second = it->second;
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)
270 ufc_integral* integral = ufc_form.create_interior_facet_integral(
id);
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);
280 if (
auto it = subdomains.find(IntegralType::interior_facet);
281 it != subdomains.end() and !interior_facet_integral_ids.empty())
283 integral_data[IntegralType::interior_facet].second = it->second;
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())
291 throw std::runtime_error(
292 "Vertex integrals not supported. Under development.");
295 return fem::Form(spaces, integral_data, coefficients, constants,
296 needs_permutation_data, mesh);
308 template <
typename T>
310 const ufc_form& ufc_form,
311 const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces,
317 const std::shared_ptr<const mesh::Mesh>& mesh =
nullptr)
320 std::vector<std::shared_ptr<const function::Function<T>>> coeff_map;
323 if (
auto it = coefficients.find(name); it != coefficients.end())
324 coeff_map.push_back(it->second);
327 throw std::runtime_error(
"Form coefficient \"" + name
328 +
"\" not provided.");
333 std::vector<std::shared_ptr<const function::Constant<T>>> const_map;
336 if (
auto it = constants.find(name); it != constants.end())
337 const_map.push_back(it->second);
340 throw std::runtime_error(
"Form constant \"" + name +
"\" not provided.");
344 return create_form(ufc_form, spaces, coeff_map, const_map, subdomains, mesh);
358 template <
typename T>
361 const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces,
367 const std::shared_ptr<const mesh::Mesh>& mesh =
nullptr)
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));
397 std::shared_ptr<function::FunctionSpace>
399 const std::string function_name,
400 std::shared_ptr<mesh::Mesh> mesh);
404 template <
typename U>
405 Eigen::Array<
typename U::scalar_type, Eigen::Dynamic, Eigen::Dynamic,
409 using T =
typename U::scalar_type;
412 const std::vector<std::shared_ptr<const function::Function<T>>> 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)
419 dofmaps[i] = coefficients[i]->function_space()->dofmap().get();
420 v.emplace_back(coefficients[i]->x()->array());
424 std::shared_ptr<const mesh::Mesh> mesh = u.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();
432 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
433 num_cells, offsets.back());
434 if (coefficients.size() > 0)
436 for (
int cell = 0; cell < num_cells; ++cell)
438 for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
440 auto dofs = dofmaps[coeff]->cell_dofs(cell);
441 const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& _v
443 for (Eigen::Index k = 0; k < dofs.size(); ++k)
444 c(cell, k + offsets[coeff]) = _v[dofs[k]];
454 template <
typename U>
455 Eigen::Array<typename U::scalar_type, Eigen::Dynamic, 1>
458 using T =
typename U::scalar_type;
460 const auto& constants = u.constants();
464 = std::accumulate(constants.begin(), constants.end(), 0,
465 [](Eigen::Index sum,
const auto& constant) {
466 return sum + constant->value.size();
470 Eigen::Array<T, Eigen::Dynamic, 1> constant_values(size);
471 Eigen::Index offset = 0;
472 for (
const auto& constant : constants)
474 const auto& value = constant->value;
475 const Eigen::Index value_size = value.size();
477 for (Eigen::Index i = 0; i < value_size; ++i)
478 constant_values[offset + i] = value[i];
480 offset += value_size;
483 return constant_values;
This class manages coordinate mappings for isoparametric cells.
Definition: CoordinateElement.h:24
A constant value which can be attached to a Form. Constants may be scalar (rank 0),...
Definition: Constant.h:19
This class represents a function in a finite element function space , given by.
Definition: Function.h:42
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:58
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
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
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
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
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
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
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
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
fem::CoordinateElement create_coordinate_map(const ufc_coordinate_mapping &ufc_cmap)
Create a CoordinateElement from ufc.
Definition: utils.cpp:228
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
IntegralType
Type of integral.
Definition: Form.h:33
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
CellType
Cell type identifier.
Definition: cell_types.h:21