10 #include <dolfinx/mesh/MeshTags.h>
51 std::vector<std::pair<
52 int, std::function<
void(T*,
const T*,
const T*,
const double*,
53 const int*,
const std::uint8_t*,
54 const std::uint32_t)>>>>& integrals,
58 for (
auto& integral_type : integrals)
60 for (
auto& integral : integral_type.second)
73 const std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
74 const std::uint8_t*,
const std::uint32_t)>&
77 return _integrals.at(
static_cast<int>(type)).at(i).tabulate;
89 const std::function<
void(T*,
const T*,
const T*,
const double*,
90 const int*,
const std::uint8_t*,
91 const std::uint32_t)>& fn)
93 std::vector<struct FormIntegrals::Integral>& integrals
94 = _integrals.at(
static_cast<int>(type));
98 for (
const auto& q : integrals)
102 throw std::runtime_error(
"Integral with ID " + std::to_string(i)
103 +
" already exists");
111 integrals.insert(integrals.begin() + pos, {fn, i, {}});
116 std::set<IntegralType>
types()
const
118 static const std::array types{IntegralType::cell,
119 IntegralType::exterior_facet,
120 IntegralType::interior_facet};
121 std::set<IntegralType> set;
122 for (
auto type : types)
123 if (!_integrals.at(
static_cast<int>(type)).empty())
133 return _integrals.at(
static_cast<int>(type)).size();
143 std::vector<int> ids;
144 for (
const auto& integral : _integrals[
static_cast<int>(type)])
145 ids.push_back(integral.id);
160 return _integrals.at(
static_cast<int>(type)).at(i).active_entities;
170 std::vector<struct Integral>& integrals
171 = _integrals.at(
static_cast<int>(type));
172 if (integrals.size() == 0)
175 std::shared_ptr<const mesh::Mesh> mesh = marker.
mesh();
177 const int tdim = topology.
dim();
179 if (type == IntegralType::exterior_facet
180 or type == IntegralType::interior_facet)
182 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
185 else if (type == IntegralType::vertex)
188 if (dim != marker.
dim())
190 throw std::runtime_error(
"Invalid MeshTags dimension:"
191 + std::to_string(marker.
dim()));
195 std::map<int, int> id_to_integral;
196 for (std::size_t i = 0; i < integrals.size(); ++i)
198 if (integrals[i].
id != -1)
200 integrals[i].active_entities.clear();
201 id_to_integral.insert({integrals[i].id, i});
206 const std::vector<int>& values = marker.
values();
207 const std::vector<std::int32_t>& tagged_entities = marker.
indices();
209 const auto entity_end
210 = std::lower_bound(tagged_entities.begin(), tagged_entities.end(),
213 if (type == IntegralType::exterior_facet)
217 if (type == IntegralType::exterior_facet)
221 std::set<std::int32_t> fwd_shared;
222 if (topology.
index_map(tdim)->num_ghosts() == 0)
225 topology.
index_map(tdim - 1)->shared_indices().begin(),
226 topology.
index_map(tdim - 1)->shared_indices().end());
229 for (
auto f = tagged_entities.begin(); f != entity_end; ++f)
231 const std::size_t i = std::distance(tagged_entities.cbegin(), f);
235 if (f_to_c->num_links(*f) == 1
236 and fwd_shared.find(*f) == fwd_shared.end())
238 if (
auto it = id_to_integral.find(values[i]);
239 it != id_to_integral.end())
241 integrals[it->second].active_entities.push_back(*f);
246 else if (type == IntegralType::interior_facet)
248 for (
auto f = tagged_entities.begin(); f != entity_end; ++f)
250 if (f_to_c->num_links(*f) == 2)
252 const std::size_t i = std::distance(tagged_entities.cbegin(), f);
253 if (
const auto it = id_to_integral.find(values[i]);
254 it != id_to_integral.end())
256 integrals[it->second].active_entities.push_back(*f);
266 for (
auto e = tagged_entities.begin(); e != entity_end; ++e)
268 const std::size_t i = std::distance(tagged_entities.cbegin(), e);
269 if (
const auto it = id_to_integral.find(values[i]);
270 it != id_to_integral.end())
272 integrals[it->second].active_entities.push_back(*e);
286 const int tdim = topology.
dim();
287 std::vector<struct Integral>& cell_integrals
288 = _integrals[
static_cast<int>(IntegralType::cell)];
291 if (cell_integrals.size() > 0 and cell_integrals[0].id == -1)
293 const int num_cells = topology.
index_map(tdim)->size_local();
294 cell_integrals[0].active_entities.resize(num_cells);
295 std::iota(cell_integrals[0].active_entities.begin(),
296 cell_integrals[0].active_entities.end(), 0);
301 std::vector<struct Integral>& exf_integrals
302 = _integrals[
static_cast<int>(IntegralType::exterior_facet)];
303 if (exf_integrals.size() > 0 and exf_integrals[0].id == -1)
306 exf_integrals[0].active_entities.clear();
312 std::set<std::int32_t> fwd_shared_facets;
315 if (topology.
index_map(tdim)->num_ghosts() == 0)
317 fwd_shared_facets.insert(
318 topology.
index_map(tdim - 1)->shared_indices().begin(),
319 topology.
index_map(tdim - 1)->shared_indices().end());
322 const int num_facets = topology.
index_map(tdim - 1)->size_local();
323 for (
int f = 0; f < num_facets; ++f)
325 if (f_to_c->num_links(f) == 1
326 and fwd_shared_facets.find(f) == fwd_shared_facets.end())
327 exf_integrals[0].active_entities.push_back(f);
333 std::vector<struct FormIntegrals::Integral>& inf_integrals
334 = _integrals[
static_cast<int>(IntegralType::interior_facet)];
335 if (!inf_integrals.empty() and inf_integrals[0].id == -1)
338 inf_integrals[0].active_entities.clear();
343 const int num_facets = topology.
index_map(tdim - 1)->size_local();
345 inf_integrals[0].active_entities.reserve(num_facets);
346 for (
int f = 0; f < num_facets; ++f)
348 if (f_to_c->num_links(f) == 2)
349 inf_integrals[0].active_entities.push_back(f);
364 std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
365 const std::uint8_t*,
const std::uint32_t)>
368 std::vector<std::int32_t> active_entities;
373 std::array<std::vector<struct Integral>, 4> _integrals;
377 bool _needs_permutation_data;