11 #include <Eigen/Dense>
12 #include <dolfinx/common/IndexMap.h>
13 #include <dolfinx/common/types.h>
14 #include <dolfinx/function/Constant.h>
15 #include <dolfinx/function/FunctionSpace.h>
16 #include <dolfinx/mesh/Geometry.h>
17 #include <dolfinx/mesh/Mesh.h>
18 #include <dolfinx/mesh/Topology.h>
22 namespace dolfinx::fem::impl
32 const mesh::Geometry& geometry,
33 const std::vector<std::int32_t>& active_cells,
34 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
35 const std::uint8_t*,
const std::uint32_t)>& fn,
36 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
38 const std::vector<T>& constant_values,
39 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info);
43 T assemble_exterior_facets(
44 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
45 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
46 const std::uint8_t*,
const std::uint32_t)>& fn,
47 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
49 const std::vector<T>& constant_values,
50 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
51 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms);
55 T assemble_interior_facets(
56 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
57 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
58 const std::uint8_t*,
const std::uint32_t)>& fn,
59 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
61 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
62 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
63 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms);
69 std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
71 const int tdim = mesh->topology().dim();
72 const std::int32_t num_cells
73 = mesh->topology().connectivity(tdim, 0)->num_nodes();
76 if (!M.all_constants_set())
77 throw std::runtime_error(
"Unset constant in Form");
78 auto constants = M.constants();
80 std::vector<T> constant_values;
81 for (
auto const& constant : constants)
84 const std::vector<T>& array = constant.second->value;
85 constant_values.insert(constant_values.end(), array.data(),
86 array.data() + array.size());
90 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
93 const FormIntegrals<T>& integrals = M.integrals();
94 const bool needs_permutation_data = integrals.needs_permutation_data();
95 if (needs_permutation_data)
96 mesh->topology_mutable().create_entity_permutations();
97 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
98 = needs_permutation_data
99 ? mesh->topology().get_cell_permutation_info()
100 : Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>(num_cells);
103 for (
int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i)
105 const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i);
106 const std::vector<std::int32_t>& active_cells
107 = integrals.integral_domains(IntegralType::cell, i);
108 value += fem::impl::assemble_cells(mesh->geometry(), active_cells, fn,
109 coeffs, constant_values, cell_info);
112 if (integrals.num_integrals(IntegralType::exterior_facet) > 0
113 or integrals.num_integrals(IntegralType::interior_facet) > 0)
116 mesh->topology_mutable().create_entities(tdim - 1);
117 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
119 const int facets_per_cell
121 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
122 = needs_permutation_data
123 ? mesh->topology().get_facet_permutations()
124 : Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>(
125 facets_per_cell, num_cells);
127 for (
int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet);
131 = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i);
132 const std::vector<std::int32_t>& active_facets
133 = integrals.integral_domains(IntegralType::exterior_facet, i);
134 value += fem::impl::assemble_exterior_facets(
135 *mesh, active_facets, fn, coeffs, constant_values, cell_info, perms);
138 const std::vector<int> c_offsets = M.coefficients().offsets();
139 for (
int i = 0; i < integrals.num_integrals(IntegralType::interior_facet);
143 = integrals.get_tabulate_tensor(IntegralType::interior_facet, i);
144 const std::vector<std::int32_t>& active_facets
145 = integrals.integral_domains(IntegralType::interior_facet, i);
146 value += fem::impl::assemble_interior_facets(
147 *mesh, active_facets, fn, coeffs, c_offsets, constant_values,
155 template <
typename T>
157 const mesh::Geometry& geometry,
158 const std::vector<std::int32_t>& active_cells,
159 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
160 const std::uint8_t*,
const std::uint32_t)>& fn,
161 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
163 const std::vector<T>& constant_values,
164 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info)
166 const int gdim = geometry.dim();
169 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
172 const int num_dofs_g = x_dofmap.num_links(0);
173 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
177 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
178 coordinate_dofs(num_dofs_g, gdim);
182 for (std::int32_t c : active_cells)
185 auto x_dofs = x_dofmap.links(c);
186 for (
int i = 0; i < num_dofs_g; ++i)
187 for (
int j = 0; j < gdim; ++j)
188 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
190 auto coeff_cell = coeffs.row(c);
191 fn(&value, coeff_cell.data(), constant_values.data(),
192 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
198 template <
typename T>
199 T assemble_exterior_facets(
200 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
201 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
202 const std::uint8_t*,
const std::uint32_t)>& fn,
203 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
205 const std::vector<T>& constant_values,
206 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
207 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
209 const int gdim = mesh.geometry().dim();
210 const int tdim = mesh.topology().dim();
213 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
216 const int num_dofs_g = x_dofmap.num_links(0);
217 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
218 = mesh.geometry().x();
221 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
222 coordinate_dofs(num_dofs_g, gdim);
224 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
226 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
231 for (std::int32_t facet : active_facets)
234 assert(f_to_c->num_links(facet) == 1);
235 const int cell = f_to_c->links(facet)[0];
238 auto facets = c_to_f->links(cell);
240 = std::find(facets.data(), facets.data() + facets.rows(), facet);
241 assert(it != (facets.data() + facets.rows()));
242 const int local_facet = std::distance(facets.data(), it);
244 auto x_dofs = x_dofmap.links(cell);
245 for (
int i = 0; i < num_dofs_g; ++i)
246 for (
int j = 0; j < gdim; ++j)
247 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
249 auto coeff_cell = coeffs.row(cell);
250 fn(&value, coeff_cell.data(), constant_values.data(),
251 coordinate_dofs.data(), &local_facet, &perms(local_facet, cell),
258 template <
typename T>
259 T assemble_interior_facets(
260 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
261 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
262 const std::uint8_t*,
const std::uint32_t)>& fn,
263 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
265 const std::vector<int>& offsets,
const std::vector<T>& constant_values,
266 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info,
267 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms)
269 const int gdim = mesh.geometry().dim();
270 const int tdim = mesh.topology().dim();
273 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
276 const int num_dofs_g = x_dofmap.num_links(0);
277 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
278 = mesh.geometry().x();
281 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
282 coordinate_dofs(2 * num_dofs_g, gdim);
283 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
284 assert(offsets.back() == coeffs.cols());
286 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
288 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
293 for (std::int32_t f : active_facets)
296 auto cells = f_to_c->links(f);
297 assert(
cells.rows() == 2);
300 std::array<int, 2> local_facet;
301 for (
int i = 0; i < 2; ++i)
303 auto facets = c_to_f->links(cells[i]);
305 = std::find(facets.data(), facets.data() + facets.rows(), f);
306 assert(it != (facets.data() + facets.rows()));
307 local_facet[i] = std::distance(facets.data(), it);
311 auto x_dofs0 = x_dofmap.links(cells[0]);
312 auto x_dofs1 = x_dofmap.links(cells[1]);
313 for (
int i = 0; i < num_dofs_g; ++i)
315 for (
int j = 0; j < gdim; ++j)
317 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
318 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
324 auto coeff_cell0 = coeffs.row(cells[0]);
325 auto coeff_cell1 = coeffs.row(cells[1]);
328 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
331 const int num_entries = offsets[i + 1] - offsets[i];
332 coeff_array.segment(2 * offsets[i], num_entries)
333 = coeff_cell0.segment(offsets[i], num_entries);
334 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
335 = coeff_cell1.segment(offsets[i], num_entries);
338 const std::array perm{perms(local_facet[0], cells[0]),
339 perms(local_facet[1], cells[1])};
340 fn(&value, coeff_array.data(), constant_values.data(),
341 coordinate_dofs.data(), local_facet.data(), perm.data(),
342 cell_info[cells[0]]);