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::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
33 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
34 const std::uint8_t*,
const std::uint32_t)>& fn,
35 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
37 const std::vector<T>& constant_values);
41 T assemble_exterior_facets(
42 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
43 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
44 const std::uint8_t*,
const std::uint32_t)>& fn,
45 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
47 const std::vector<T>& constant_values);
51 T assemble_interior_facets(
52 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
53 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
54 const std::uint8_t*,
const std::uint32_t)>& fn,
55 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
57 const std::vector<int>& offsets,
const std::vector<T>& constant_values);
63 std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
67 if (!M.all_constants_set())
68 throw std::runtime_error(
"Unset constant in Form");
69 auto constants = M.constants();
71 std::vector<T> constant_values;
72 for (
auto const& constant : constants)
75 const std::vector<T>& array = constant.second->value;
76 constant_values.insert(constant_values.end(), array.data(),
77 array.data() + array.size());
81 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> coeffs
84 const FormIntegrals<T>& integrals = M.integrals();
86 for (
int i = 0; i < integrals.num_integrals(IntegralType::cell); ++i)
88 const auto& fn = integrals.get_tabulate_tensor(IntegralType::cell, i);
89 const std::vector<std::int32_t>& active_cells
90 = integrals.integral_domains(IntegralType::cell, i);
91 value += fem::impl::assemble_cells(*mesh, active_cells, fn, coeffs,
95 for (
int i = 0; i < integrals.num_integrals(IntegralType::exterior_facet);
99 = integrals.get_tabulate_tensor(IntegralType::exterior_facet, i);
100 const std::vector<std::int32_t>& active_facets
101 = integrals.integral_domains(IntegralType::exterior_facet, i);
102 value += fem::impl::assemble_exterior_facets(*mesh, active_facets, fn,
103 coeffs, constant_values);
106 const std::vector<int> c_offsets = M.coefficients().offsets();
107 for (
int i = 0; i < integrals.num_integrals(IntegralType::interior_facet);
111 = integrals.get_tabulate_tensor(IntegralType::interior_facet, i);
112 const std::vector<std::int32_t>& active_facets
113 = integrals.integral_domains(IntegralType::interior_facet, i);
114 value += fem::impl::assemble_interior_facets(
115 *mesh, active_facets, fn, coeffs, c_offsets, constant_values);
121 template <
typename T>
123 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_cells,
124 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
125 const std::uint8_t*,
const std::uint32_t)>& fn,
126 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
128 const std::vector<T>& constant_values)
130 const int gdim = mesh.geometry().dim();
131 const int tdim = mesh.topology().dim();
132 mesh.topology_mutable().create_entities(tdim);
133 mesh.topology_mutable().create_entity_permutations();
136 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
139 const int num_dofs_g = x_dofmap.num_links(0);
140 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
141 = mesh.geometry().x();
144 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
145 coordinate_dofs(num_dofs_g, gdim);
147 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
148 = mesh.topology().get_cell_permutation_info();
152 for (std::int32_t c : active_cells)
155 auto x_dofs = x_dofmap.links(c);
156 for (
int i = 0; i < num_dofs_g; ++i)
157 for (
int j = 0; j < gdim; ++j)
158 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
160 auto coeff_cell = coeffs.row(c);
161 fn(&value, coeff_cell.data(), constant_values.data(),
162 coordinate_dofs.data(),
nullptr,
nullptr, cell_info[c]);
168 template <
typename T>
169 T assemble_exterior_facets(
170 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
171 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
172 const std::uint8_t*,
const std::uint32_t)>& fn,
173 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
175 const std::vector<T>& constant_values)
177 const int gdim = mesh.geometry().dim();
178 const int tdim = mesh.topology().dim();
181 mesh.topology_mutable().create_entities(tdim - 1);
182 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
183 mesh.topology_mutable().create_entity_permutations();
186 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
189 const int num_dofs_g = x_dofmap.num_links(0);
190 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
191 = mesh.geometry().x();
194 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
195 coordinate_dofs(num_dofs_g, gdim);
197 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
198 = mesh.topology().get_facet_permutations();
199 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
200 = mesh.topology().get_cell_permutation_info();
202 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
204 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
209 for (std::int32_t facet : active_facets)
212 assert(f_to_c->num_links(facet) == 1);
213 const int cell = f_to_c->links(facet)[0];
216 auto facets = c_to_f->links(cell);
218 = std::find(facets.data(), facets.data() + facets.rows(), facet);
219 assert(it != (facets.data() + facets.rows()));
220 const int local_facet = std::distance(facets.data(), it);
222 auto x_dofs = x_dofmap.links(cell);
223 for (
int i = 0; i < num_dofs_g; ++i)
224 for (
int j = 0; j < gdim; ++j)
225 coordinate_dofs(i, j) = x_g(x_dofs[i], j);
227 auto coeff_cell = coeffs.row(cell);
228 const std::uint8_t perm = perms(local_facet, cell);
229 fn(&value, coeff_cell.data(), constant_values.data(),
230 coordinate_dofs.data(), &local_facet, &perm, cell_info[cell]);
236 template <
typename T>
237 T assemble_interior_facets(
238 const mesh::Mesh& mesh,
const std::vector<std::int32_t>& active_facets,
239 const std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
240 const std::uint8_t*,
const std::uint32_t)>& fn,
241 const Eigen::Array<T, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>&
243 const std::vector<int>& offsets,
const std::vector<T>& constant_values)
245 const int gdim = mesh.geometry().dim();
246 const int tdim = mesh.topology().dim();
249 mesh.topology_mutable().create_entities(tdim - 1);
250 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
251 mesh.topology_mutable().create_entity_permutations();
254 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
257 const int num_dofs_g = x_dofmap.num_links(0);
258 const Eigen::Array<double, Eigen::Dynamic, 3, Eigen::RowMajor>& x_g
259 = mesh.geometry().x();
262 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
263 coordinate_dofs(2 * num_dofs_g, gdim);
264 Eigen::Array<T, Eigen::Dynamic, 1> coeff_array(2 * offsets.back());
265 assert(offsets.back() == coeffs.cols());
267 const Eigen::Array<std::uint8_t, Eigen::Dynamic, Eigen::Dynamic>& perms
268 = mesh.topology().get_facet_permutations();
269 const Eigen::Array<std::uint32_t, Eigen::Dynamic, 1>& cell_info
270 = mesh.topology().get_cell_permutation_info();
272 auto f_to_c = mesh.topology().connectivity(tdim - 1, tdim);
274 auto c_to_f = mesh.topology().connectivity(tdim, tdim - 1);
279 for (std::int32_t f : active_facets)
282 auto cells = f_to_c->links(f);
283 assert(cells.rows() == 2);
286 std::array<int, 2> local_facet;
287 for (
int i = 0; i < 2; ++i)
289 auto facets = c_to_f->links(cells[i]);
291 = std::find(facets.data(), facets.data() + facets.rows(), f);
292 assert(it != (facets.data() + facets.rows()));
293 local_facet[i] = std::distance(facets.data(), it);
296 const std::array perm{perms(local_facet[0], cells[0]),
297 perms(local_facet[1], cells[1])};
300 auto x_dofs0 = x_dofmap.links(cells[0]);
301 auto x_dofs1 = x_dofmap.links(cells[1]);
302 for (
int i = 0; i < num_dofs_g; ++i)
304 for (
int j = 0; j < gdim; ++j)
306 coordinate_dofs(i, j) = x_g(x_dofs0[i], j);
307 coordinate_dofs(i + num_dofs_g, j) = x_g(x_dofs1[i], j);
313 auto coeff_cell0 = coeffs.row(cells[0]);
314 auto coeff_cell1 = coeffs.row(cells[1]);
317 for (std::size_t i = 0; i < offsets.size() - 1; ++i)
320 const int num_entries = offsets[i + 1] - offsets[i];
321 coeff_array.segment(2 * offsets[i], num_entries)
322 = coeff_cell0.segment(offsets[i], num_entries);
323 coeff_array.segment(offsets[i + 1] + offsets[i], num_entries)
324 = coeff_cell1.segment(offsets[i], num_entries);
326 fn(&value, coeff_array.data(), constant_values.data(),
327 coordinate_dofs.data(), local_facet.data(), perm.data(),
328 cell_info[cells[0]]);