9 #include <dolfinx/function/FunctionSpace.h>
10 #include <dolfinx/mesh/Mesh.h>
11 #include <dolfinx/mesh/MeshTags.h>
81 Form(
const std::vector<std::shared_ptr<const function::FunctionSpace>>&
86 std::vector<std::pair<
87 int, std::function<
void(
88 T*,
const T*,
const T*,
const double*,
const int*,
89 const std::uint8_t*,
const std::uint32_t)>>>,
96 const std::shared_ptr<const mesh::Mesh>&
mesh =
nullptr)
105 if (_mesh != V->mesh())
106 throw std::runtime_error(
"Incompatible mesh");
108 throw std::runtime_error(
"No mesh could be associated with the Form.");
111 for (
auto& integral_type : integrals)
115 auto it = _integrals.emplace(
116 type, std::map<
int, std::pair<kern, std::vector<std::int32_t>>>());
119 for (
auto& integral : integral_type.second.first)
120 it.first->second.insert({integral.first, {integral.second, {}}});
124 if (integral_type.second.second)
126 assert(_mesh == integral_type.second.second->mesh());
127 set_domains(type, *integral_type.second.second);
133 set_default_domains(*_mesh);
148 int rank()
const {
return _function_spaces.size(); }
152 std::shared_ptr<const mesh::Mesh>
mesh()
const {
return _mesh; }
156 const std::vector<std::shared_ptr<const function::FunctionSpace>>&
159 return _function_spaces;
167 const std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
168 const std::uint8_t*,
const std::uint32_t)>&
171 auto it0 = _integrals.find(type);
172 if (it0 == _integrals.end())
173 throw std::runtime_error(
"No kernels for requested type.");
174 auto it1 = it0->second.find(i);
175 if (it1 == it0->second.end())
176 throw std::runtime_error(
"No kernel for requested domain index.");
178 return it1->second.first;
185 std::set<IntegralType> set;
186 for (
auto& type : _integrals)
187 set.insert(type.first);
196 if (
auto it = _integrals.find(type); it == _integrals.end())
199 return it->second.size();
210 std::vector<int> ids;
211 if (
auto it = _integrals.find(type); it != _integrals.end())
213 for (
auto&
kernel : it->second)
214 ids.push_back(
kernel.first);
227 auto it0 = _integrals.find(type);
228 if (it0 == _integrals.end())
229 throw std::runtime_error(
"No mesh entities for requested type.");
230 auto it1 = it0->second.find(i);
231 if (it1 == it0->second.end())
232 throw std::runtime_error(
"No mesh entities for requested domain index.");
233 return it1->second.second;
237 const std::vector<std::shared_ptr<const function::Function<T>>>
240 return _coefficients;
253 std::vector<int> n{0};
254 for (
const auto& c : _coefficients)
257 throw std::runtime_error(
"Not all form coefficients have been set.");
258 n.push_back(n.back() + c->function_space()->element()->space_dimension());
264 const std::vector<std::shared_ptr<const function::Constant<T>>>&
281 auto it0 = _integrals.find(type);
282 assert(it0 != _integrals.end());
284 std::shared_ptr<const mesh::Mesh>
mesh = marker.
mesh();
286 const int tdim = topology.
dim();
288 if (type == IntegralType::exterior_facet
289 or type == IntegralType::interior_facet)
292 mesh->topology_mutable().create_connectivity(dim, tdim);
294 else if (type == IntegralType::vertex)
297 if (dim != marker.
dim())
299 throw std::runtime_error(
"Invalid MeshTags dimension:"
300 + std::to_string(marker.
dim()));
304 std::map<int, std::pair<kern, std::vector<std::int32_t>>>& integrals
308 const std::vector<int>& values = marker.
values();
309 const std::vector<std::int32_t>& tagged_entities = marker.
indices();
311 const auto entity_end
312 = std::lower_bound(tagged_entities.begin(), tagged_entities.end(),
319 if (type == IntegralType::exterior_facet)
324 std::set<std::int32_t> fwd_shared;
325 if (topology.
index_map(tdim)->num_ghosts() == 0)
328 topology.
index_map(tdim - 1)->shared_indices().begin(),
329 topology.
index_map(tdim - 1)->shared_indices().end());
332 for (
auto f = tagged_entities.begin(); f != entity_end; ++f)
336 if (f_to_c->num_links(*f) == 1
337 and fwd_shared.find(*f) == fwd_shared.end())
339 const std::size_t i = std::distance(tagged_entities.cbegin(), f);
340 if (
auto it = integrals.find(values[i]); it != integrals.end())
341 it->second.second.push_back(*f);
345 else if (type == IntegralType::interior_facet)
347 for (
auto f = tagged_entities.begin(); f != entity_end; ++f)
349 if (f_to_c->num_links(*f) == 2)
351 const std::size_t i = std::distance(tagged_entities.cbegin(), f);
352 if (
auto it = integrals.find(values[i]); it != integrals.end())
353 it->second.second.push_back(*f);
362 for (
auto e = tagged_entities.begin(); e != entity_end; ++e)
364 const std::size_t i = std::distance(tagged_entities.cbegin(), e);
365 if (
auto it = integrals.find(values[i]); it != integrals.end())
366 it->second.second.push_back(*e);
379 const int tdim = topology.
dim();
382 if (
auto kernels = _integrals.find(IntegralType::cell);
383 kernels != _integrals.end())
385 if (
auto it = kernels->second.find(-1); it != kernels->second.end())
387 std::vector<std::int32_t>& active_entities = it->second.second;
388 const int num_cells = topology.
index_map(tdim)->size_local();
389 active_entities.resize(num_cells);
390 std::iota(active_entities.begin(), active_entities.end(), 0);
396 if (
auto kernels = _integrals.find(IntegralType::exterior_facet);
397 kernels != _integrals.end())
399 if (
auto it = kernels->second.find(-1); it != kernels->second.end())
401 std::vector<std::int32_t>& active_entities = it->second.second;
402 active_entities.clear();
405 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
408 std::set<std::int32_t> fwd_shared_facets;
411 if (topology.
index_map(tdim)->num_ghosts() == 0)
413 fwd_shared_facets.insert(
414 topology.
index_map(tdim - 1)->shared_indices().begin(),
415 topology.
index_map(tdim - 1)->shared_indices().end());
418 const int num_facets = topology.
index_map(tdim - 1)->size_local();
419 for (
int f = 0; f < num_facets; ++f)
421 if (f_to_c->num_links(f) == 1
422 and fwd_shared_facets.find(f) == fwd_shared_facets.end())
424 active_entities.push_back(f);
432 if (
auto kernels = _integrals.find(IntegralType::interior_facet);
433 kernels != _integrals.end())
435 if (
auto it = kernels->second.find(-1); it != kernels->second.end())
437 std::vector<std::int32_t>& active_entities = it->second.second;
440 mesh.topology_mutable().create_connectivity(tdim - 1, tdim);
442 const int num_facets = topology.
index_map(tdim - 1)->size_local();
444 active_entities.clear();
445 active_entities.reserve(num_facets);
446 for (
int f = 0; f < num_facets; ++f)
448 if (f_to_c->num_links(f) == 2)
449 active_entities.push_back(f);
456 std::vector<std::shared_ptr<const function::FunctionSpace>> _function_spaces;
459 std::vector<std::shared_ptr<const function::Function<T>>> _coefficients;
462 std::vector<std::shared_ptr<const function::Constant<T>>> _constants;
465 std::shared_ptr<const mesh::Mesh> _mesh;
468 = std::function<void(T*,
const T*,
const T*,
const double*,
const int*,
469 const std::uint8_t*,
const std::uint32_t)>;
471 std::map<int, std::pair<kern, std::vector<std::int32_t>>>>
475 bool _needs_permutation_data;
A constant value which can be attached to a Form. Constants may be scalar (rank 0),...
Definition: Constant.h:19
This class represents a function in a finite element function space , given by.
Definition: Function.h:42
A Mesh consists of a set of connected and numbered mesh topological entities, and geometry data.
Definition: Mesh.h:47
Topology stores the topology of a mesh, consisting of mesh entities and connectivity (incidence relat...
Definition: Topology.h:58
std::shared_ptr< const graph::AdjacencyList< std::int32_t > > connectivity(int d0, int d1) const
Return connectivity from entities of dimension d0 to entities of dimension d1.
Definition: Topology.cpp:255
std::shared_ptr< const common::IndexMap > index_map(int dim) const
Get the IndexMap that described the parallel distribution of the mesh entities.
Definition: Topology.cpp:162
int dim() const
Return topological dimension.
Definition: Topology.cpp:153
IntegralType
Type of integral.
Definition: Form.h:33