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>
25 struct ufc_coordinate_mapping;
26 struct ufc_function_space;
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)
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)
73 spaces(i, j) = {a(i, j)->function_space(0), a(i, j)->function_space(1)};
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);
102 throw std::runtime_error(
103 "Cannot create sparsity pattern. Form is not a bilinear form");
109 std::shared_ptr mesh = a.
mesh();
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())
117 const int tdim = mesh->topology().dim();
118 mesh->topology_mutable().create_entities(tdim - 1);
119 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
130 const std::array<const DofMap*, 2>& dofmaps,
131 const std::set<IntegralType>& integrals);
136 const std::vector<int>& parent_map
143 DofMap
create_dofmap(MPI_Comm comm,
const ufc_dofmap& dofmap,
144 mesh::Topology& topology);
147 template <
typename T>
149 std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
153 std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
155 const char** names = ufc_form.coefficient_name_map();
156 for (
int i = 0; i < ufc_form.num_coefficients; ++i)
158 coeffs.emplace_back(ufc_form.original_coefficient_position(i), names[i],
165 template <
typename T>
167 std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
171 std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
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);
182 template <
typename T>
184 const ufc_form& ufc_form,
185 const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
187 assert(ufc_form.rank == (
int)spaces.size());
190 for (std::size_t i = 0; i < spaces.size(); ++i)
192 assert(spaces[i]->element());
193 std::unique_ptr<ufc_finite_element, decltype(free)*> ufc_element(
194 ufc_form.create_finite_element(i), free);
196 if (std::string(ufc_element->signature)
197 != spaces[i]->element()->signature())
199 throw std::runtime_error(
200 "Cannot create form. Wrong type of function space for argument.");
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)
211 ufc_integral* cell_integral = ufc_form.create_cell_integral(
id);
212 assert(cell_integral);
214 cell_integral->tabulate_tensor);
215 std::free(cell_integral);
220 if (ufc_form.num_exterior_facet_integrals > 0
221 or ufc_form.num_interior_facet_integrals > 0)
225 auto mesh = spaces[0]->mesh();
226 const int tdim = mesh->topology().dim();
227 spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
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)
236 ufc_integral* exterior_facet_integral
237 = ufc_form.create_exterior_facet_integral(
id);
238 assert(exterior_facet_integral);
240 exterior_facet_integral->tabulate_tensor);
241 std::free(exterior_facet_integral);
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)
249 ufc_integral* interior_facet_integral
250 = ufc_form.create_interior_facet_integral(
id);
251 assert(interior_facet_integral);
253 interior_facet_integral->tabulate_tensor);
254 std::free(interior_facet_integral);
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)
262 throw std::runtime_error(
263 "Vertex integrals not supported. Under development.");
269 fem::get_constants_from_ufc_form<T>(ufc_form));
278 template <
typename T>
281 const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
283 ufc_form* form = fptr();
284 auto L = std::make_shared<fem::Form<T>>(
285 dolfinx::fem::create_form<T>(*form, spaces));
311 std::shared_ptr<function::FunctionSpace>
313 const std::string function_name,
314 std::shared_ptr<mesh::Mesh> mesh);
318 template <
typename T>
319 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
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)
329 dofmaps[i] = coefficients.
get(i)->function_space()->dofmap().get();
330 v.emplace_back(coefficients.
get(i)->x()->array());
334 std::shared_ptr<const mesh::Mesh> mesh = form.
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();
342 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
343 num_cells, offsets.back());
344 if (coefficients.
size() > 0)
346 for (
int cell = 0; cell < num_cells; ++cell)
348 for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
350 auto dofs = dofmaps[coeff]->cell_dofs(cell);
351 const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& _v
353 for (Eigen::Index k = 0; k < dofs.size(); ++k)
354 c(cell, k + offsets[coeff]) = _v[dofs[k]];
364 template <
typename T>
367 std::vector<T> constant_values;
370 const std::vector<T>& array = constant.second->value;
371 constant_values.insert(constant_values.end(), array.begin(), array.end());
374 return Eigen::Map<const Eigen::Array<T, Eigen::Dynamic, 1>>(
375 constant_values.data(), constant_values.size(), 1);