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>
58 std::vector<std::array<std::shared_ptr<const function::FunctionSpace>, 2>>>
60 const Eigen::Ref<
const Eigen::Array<
const fem::Form<T>*, Eigen::Dynamic,
61 Eigen::Dynamic, Eigen::RowMajor>>& a)
63 std::vector<std::vector<
64 std::array<std::shared_ptr<const function::FunctionSpace>, 2>>>
67 std::array<std::shared_ptr<const function::FunctionSpace>, 2>>(
69 for (
int i = 0; i < a.rows(); ++i)
70 for (
int j = 0; j < a.cols(); ++j)
72 spaces[i][j] = {a(i, j)->function_space(0), a(i, j)->function_space(1)};
86 throw std::runtime_error(
87 "Cannot create sparsity pattern. Form is not a bilinear form");
93 std::shared_ptr mesh = a.
mesh();
96 const std::set<IntegralType> types = a.
integrals().types();
97 if (types.find(IntegralType::interior_facet) != types.end()
98 or types.find(IntegralType::exterior_facet) != types.end())
101 const int tdim = mesh->topology().dim();
102 mesh->topology_mutable().create_entities(tdim - 1);
103 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
114 const std::array<const DofMap*, 2>& dofmaps,
115 const std::set<IntegralType>& integrals);
120 const std::vector<int>& parent_map
127 DofMap
create_dofmap(MPI_Comm comm,
const ufc_dofmap& dofmap,
128 mesh::Topology& topology);
131 template <
typename T>
133 std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
137 std::tuple<int, std::string, std::shared_ptr<function::Function<T>>>>
139 const char** names = ufc_form.coefficient_name_map();
140 for (
int i = 0; i < ufc_form.num_coefficients; ++i)
142 coeffs.emplace_back(ufc_form.original_coefficient_position(i), names[i],
149 template <
typename T>
151 std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
155 std::pair<std::string, std::shared_ptr<const function::Constant<T>>>>
157 const char** names = ufc_form.constant_name_map();
158 for (
int i = 0; i < ufc_form.num_constants; ++i)
159 constants.emplace_back(names[i],
nullptr);
166 template <
typename T>
168 const ufc_form& ufc_form,
169 const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
171 assert(ufc_form.rank == (
int)spaces.size());
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.");
189 bool needs_permutation_data =
false;
191 std::vector<std::pair<
192 int, std::function<void(T*,
const T*,
const T*,
const double*,
193 const int*,
const std::uint8_t*,
194 const std::uint32_t)>>>>
197 std::vector<int> cell_integral_ids(ufc_form.num_cell_integrals);
198 ufc_form.get_cell_integral_ids(cell_integral_ids.data());
199 for (
int id : cell_integral_ids)
201 ufc_integral* integral = ufc_form.create_cell_integral(
id);
203 if (integral->needs_permutation_data)
204 needs_permutation_data =
true;
205 integral_data[IntegralType::cell].emplace_back(
id,
206 integral->tabulate_tensor);
212 if (ufc_form.num_exterior_facet_integrals > 0
213 or ufc_form.num_interior_facet_integrals > 0)
217 auto mesh = spaces[0]->mesh();
218 const int tdim = mesh->topology().dim();
219 spaces[0]->mesh()->topology_mutable().create_entities(tdim - 1);
223 std::vector<int> exterior_facet_integral_ids(
224 ufc_form.num_exterior_facet_integrals);
225 ufc_form.get_exterior_facet_integral_ids(exterior_facet_integral_ids.data());
226 for (
int id : exterior_facet_integral_ids)
228 ufc_integral* integral = ufc_form.create_exterior_facet_integral(
id);
230 if (integral->needs_permutation_data)
231 needs_permutation_data =
true;
232 integral_data[IntegralType::exterior_facet].emplace_back(
233 id, integral->tabulate_tensor);
237 std::vector<int> interior_facet_integral_ids(
238 ufc_form.num_interior_facet_integrals);
239 ufc_form.get_interior_facet_integral_ids(interior_facet_integral_ids.data());
240 for (
int id : interior_facet_integral_ids)
242 ufc_integral* integral = ufc_form.create_interior_facet_integral(
id);
244 if (integral->needs_permutation_data)
245 needs_permutation_data =
true;
246 integral_data[IntegralType::interior_facet].emplace_back(
247 id, integral->tabulate_tensor);
252 std::vector<int> vertex_integral_ids(ufc_form.num_vertex_integrals);
253 ufc_form.get_vertex_integral_ids(vertex_integral_ids.data());
254 if (vertex_integral_ids.size() > 0)
256 throw std::runtime_error(
257 "Vertex integrals not supported. Under development.");
263 fem::get_constants_from_ufc_form<T>(ufc_form));
272 template <
typename T>
275 const std::vector<std::shared_ptr<const function::FunctionSpace>>& spaces)
277 ufc_form* form = fptr();
278 auto L = std::make_shared<fem::Form<T>>(
279 dolfinx::fem::create_form<T>(*form, spaces));
305 std::shared_ptr<function::FunctionSpace>
307 const std::string function_name,
308 std::shared_ptr<mesh::Mesh> mesh);
312 template <
typename T>
313 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
318 const std::vector<int>& offsets = coefficients.
offsets();
319 std::vector<const fem::DofMap*> dofmaps(coefficients.
size());
320 std::vector<Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>> v;
321 for (
int i = 0; i < coefficients.
size(); ++i)
323 dofmaps[i] = coefficients.
get(i)->function_space()->dofmap().get();
324 v.emplace_back(coefficients.
get(i)->x()->array());
328 std::shared_ptr<const mesh::Mesh> mesh = form.
mesh();
330 const int tdim = mesh->topology().dim();
331 const std::int32_t num_cells
332 = mesh->topology().index_map(tdim)->size_local()
333 + mesh->topology().index_map(tdim)->num_ghosts();
336 Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> c(
337 num_cells, offsets.back());
338 if (coefficients.
size() > 0)
340 for (
int cell = 0; cell < num_cells; ++cell)
342 for (std::size_t coeff = 0; coeff < dofmaps.size(); ++coeff)
344 auto dofs = dofmaps[coeff]->cell_dofs(cell);
345 const Eigen::Ref<const Eigen::Matrix<T, Eigen::Dynamic, 1>>& _v
347 for (Eigen::Index k = 0; k < dofs.size(); ++k)
348 c(cell, k + offsets[coeff]) = _v[dofs[k]];
358 template <
typename T>
361 std::vector<T> constant_values;
364 const std::vector<T>& array = constant.second->value;
365 constant_values.insert(constant_values.end(), array.begin(), array.end());
368 return Eigen::Map<const Eigen::Array<T, Eigen::Dynamic, 1>>(
369 constant_values.data(), constant_values.size(), 1);