10 #include <dolfinx/mesh/MeshTags.h>
49 const std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
50 const std::uint8_t*,
const std::uint32_t)>&
53 return _integrals.at(
static_cast<int>(type)).at(i).tabulate;
63 std::function<
void(T*,
const T*,
const T*,
const double*,
const int*,
64 const std::uint8_t*,
const std::uint32_t)>
67 std::vector<struct FormIntegrals::Integral>& integrals
68 = _integrals.at(
static_cast<int>(type));
72 for (
const auto& q : integrals)
76 throw std::runtime_error(
"Integral with ID " + std::to_string(i)
85 integrals.insert(integrals.begin() + pos, {fn, i, {}});
90 std::set<IntegralType>
types()
const
92 static const std::array types{IntegralType::cell,
93 IntegralType::exterior_facet,
94 IntegralType::interior_facet};
95 std::set<IntegralType> set;
96 for (
auto type : types)
97 if (!_integrals.at(
static_cast<int>(type)).empty())
107 return _integrals.at(
static_cast<int>(type)).size();
117 std::vector<int> ids;
118 for (
const auto& integral : _integrals[
static_cast<int>(type)])
119 ids.push_back(integral.id);
134 return _integrals.at(
static_cast<int>(type)).at(i).active_entities;
144 std::vector<struct Integral>& integrals
145 = _integrals.at(
static_cast<int>(type));
146 if (integrals.size() == 0)
149 std::shared_ptr<const mesh::Mesh> mesh = marker.
mesh();
151 const int tdim = topology.
dim();
153 if (type == IntegralType::exterior_facet
154 or type == IntegralType::interior_facet)
156 mesh->topology_mutable().create_connectivity(tdim - 1, tdim);
159 else if (type == IntegralType::vertex)
162 if (dim != marker.
dim())
164 throw std::runtime_error(
"Invalid MeshTags dimension:"
165 + std::to_string(marker.
dim()));
169 std::map<int, int> id_to_integral;
170 for (std::size_t i = 0; i < integrals.size(); ++i)
172 if (integrals[i].
id != -1)
174 integrals[i].active_entities.clear();
175 id_to_integral.insert({integrals[i].id, i});
180 const std::vector<int>& values = marker.
values();
181 const std::vector<std::int32_t>& tagged_entities = marker.
indices();
183 const auto entity_end
184 = std::lower_bound(tagged_entities.begin(), tagged_entities.end(),
187 if (type == IntegralType::exterior_facet)
191 if (type == IntegralType::exterior_facet)
195 std::set<std::int32_t> fwd_shared;
196 if (topology.
index_map(tdim)->num_ghosts() == 0)
199 topology.
index_map(tdim - 1)->shared_indices().begin(),
200 topology.
index_map(tdim - 1)->shared_indices().end());
203 for (
auto f = tagged_entities.begin(); f != entity_end; ++f)
205 const std::size_t i = std::distance(tagged_entities.cbegin(), f);
209 if (f_to_c->num_links(*f) == 1
210 and fwd_shared.find(*f) == fwd_shared.end())
212 if (
auto it = id_to_integral.find(values[i]);
213 it != id_to_integral.end())
215 integrals[it->second].active_entities.push_back(*f);
220 else if (type == IntegralType::interior_facet)
222 for (
auto f = tagged_entities.begin(); f != entity_end; ++f)
224 if (f_to_c->num_links(*f) == 2)
226 const std::size_t i = std::distance(tagged_entities.cbegin(), f);
227 if (
const auto it = id_to_integral.find(values[i]);
228 it != id_to_integral.end())
230 integrals[it->second].active_entities.push_back(*f);
240 for (
auto e = tagged_entities.begin(); e != entity_end; ++e)
242 const std::size_t i = std::distance(tagged_entities.cbegin(), e);
243 if (
const auto it = id_to_integral.find(values[i]);
244 it != id_to_integral.end())
246 integrals[it->second].active_entities.push_back(*e);
260 const int tdim = topology.
dim();
261 std::vector<struct Integral>& cell_integrals
262 = _integrals[
static_cast<int>(IntegralType::cell)];
265 if (cell_integrals.size() > 0 and cell_integrals[0].id == -1)
267 const int num_cells = topology.
index_map(tdim)->size_local();
268 cell_integrals[0].active_entities.resize(num_cells);
269 std::iota(cell_integrals[0].active_entities.begin(),
270 cell_integrals[0].active_entities.end(), 0);
275 std::vector<struct Integral>& exf_integrals
276 = _integrals[
static_cast<int>(IntegralType::exterior_facet)];
277 if (exf_integrals.size() > 0 and exf_integrals[0].id == -1)
280 exf_integrals[0].active_entities.clear();
286 std::set<std::int32_t> fwd_shared_facets;
289 if (topology.
index_map(tdim)->num_ghosts() == 0)
291 fwd_shared_facets.insert(
292 topology.
index_map(tdim - 1)->shared_indices().begin(),
293 topology.
index_map(tdim - 1)->shared_indices().end());
296 const int num_facets = topology.
index_map(tdim - 1)->size_local();
297 for (
int f = 0; f < num_facets; ++f)
299 if (f_to_c->num_links(f) == 1
300 and fwd_shared_facets.find(f) == fwd_shared_facets.end())
301 exf_integrals[0].active_entities.push_back(f);
307 std::vector<struct FormIntegrals::Integral>& inf_integrals
308 = _integrals[
static_cast<int>(IntegralType::interior_facet)];
309 if (!inf_integrals.empty() and inf_integrals[0].id == -1)
312 inf_integrals[0].active_entities.clear();
317 const int num_facets = topology.
index_map(tdim - 1)->size_local();
319 inf_integrals[0].active_entities.reserve(num_facets);
320 for (
int f = 0; f < num_facets; ++f)
322 if (f_to_c->num_links(f) == 2)
323 inf_integrals[0].active_entities.push_back(f);
333 std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
334 const std::uint8_t*,
const std::uint32_t)>
337 std::vector<std::int32_t> active_entities;
342 std::array<std::vector<struct Integral>, 4> _integrals;