10 #include <dolfinx/common/IndexMap.h>
11 #include <dolfinx/fem/DofMap.h>
12 #include <dolfinx/fem/FunctionSpace.h>
13 #include <dolfinx/la/SparsityPattern.h>
14 #include <dolfinx/mesh/Mesh.h>
62 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
63 const std::int32_t*,
const T*)>& mat_set,
67 using namespace dolfinx;
70 const std::function<
int(std::int32_t,
const std::int32_t*, std::int32_t,
71 const std::int32_t*,
const T*)>& mat_set,
75 std::shared_ptr<const mesh::Mesh> mesh = V0.
mesh();
80 if (mesh != V1.
mesh())
82 throw std::runtime_error(
"Compute discrete gradient operator. Function "
83 "spaces do not share the same mesh");
87 mesh->topology_mutable().create_entities(1);
88 std::int64_t num_edges_global = mesh->topology().index_map(1)->size_global();
89 const std::int64_t V0dim
90 = V0.
dofmap()->index_map->size_global() * V0.
dofmap()->index_map_bs();
91 if (V0dim != num_edges_global)
93 throw std::runtime_error(
94 "Cannot compute discrete gradient operator. Function "
95 "spaces is not a lowest-order edge space");
99 const std::int64_t num_vertices_global
100 = mesh->topology().index_map(0)->size_global();
101 const std::int64_t V1dim
102 = V1.
dofmap()->index_map->size_global() * V1.
dofmap()->index_map_bs();
103 if (V1dim != num_vertices_global)
105 throw std::runtime_error(
106 "Cannot compute discrete gradient operator. Function "
107 "space is not a linear nodal function space");
111 std::shared_ptr<const dolfinx::fem::ElementDofLayout> layout0
112 = V0.
dofmap()->element_dof_layout;
113 std::shared_ptr<const dolfinx::fem::ElementDofLayout> layout1
114 = V1.
dofmap()->element_dof_layout;
117 std::array<std::shared_ptr<const common::IndexMap>, 2> index_maps
119 std::array<int, 2> block_sizes
120 = {V0.
dofmap()->index_map_bs(), V1.
dofmap()->index_map_bs()};
121 std::vector<std::array<std::int64_t, 2>> local_range
122 = {index_maps[0]->local_range(), index_maps[1]->local_range()};
123 assert(block_sizes[0] == block_sizes[1]);
126 const int tdim = mesh->topology().dim();
127 mesh->topology_mutable().create_connectivity(1, 0);
128 auto e_to_v = mesh->topology().connectivity(1, 0);
129 mesh->topology_mutable().create_connectivity(tdim, 1);
130 auto c_to_e = mesh->topology().connectivity(tdim, 1);
131 mesh->topology_mutable().create_connectivity(1, tdim);
132 auto e_to_c = mesh->topology().connectivity(1, tdim);
133 mesh->topology_mutable().create_connectivity(tdim, 0);
134 auto c_to_v = mesh->topology().connectivity(tdim, 0);
137 const std::int32_t num_edges = mesh->topology().index_map(1)->size_local()
138 + mesh->topology().index_map(1)->num_ghosts();
139 const std::shared_ptr<const fem::DofMap> dofmap0 = V0.
dofmap();
142 const int num_edges_per_cell
144 std::map<std::int32_t, std::vector<std::int32_t>> local_edge_dofs;
145 for (std::int32_t i = 0; i < num_edges_per_cell; ++i)
146 local_edge_dofs[i] = layout0->entity_dofs(1, i);
148 const int num_vertices_per_cell
150 std::map<std::int32_t, std::vector<std::int32_t>> local_vertex_dofs;
151 for (std::int32_t i = 0; i < num_vertices_per_cell; ++i)
152 local_vertex_dofs[i] = layout1->entity_dofs(0, i);
155 const std::shared_ptr<const fem::DofMap> dofmap1 = V1.
dofmap();
157 const std::vector<std::int64_t>& global_indices
158 = mesh->topology().index_map(0)->global_indices();
160 for (std::int32_t e = 0; e < num_edges; ++e)
163 tcb::span<const std::int32_t> cells = e_to_c->links(e);
164 assert(cells.size() > 0);
165 const std::int32_t cell = cells[0];
166 tcb::span<const std::int32_t> edges = c_to_e->links(cell);
167 const auto it = std::find(edges.begin(), edges.end(), e);
168 assert(it != edges.end());
169 const int local_edge = std::distance(edges.begin(), it);
172 tcb::span<const std::int32_t> dofs0 = dofmap0->cell_dofs(cell);
173 std::vector<std::int32_t>& local_dofs = local_edge_dofs[local_edge];
174 assert(local_dofs.size() == 1);
175 const std::int32_t row = dofs0[local_dofs[0]];
177 tcb::span<const std::int32_t> vertices = e_to_v->links(e);
178 assert(vertices.size() == 2);
179 tcb::span<const std::int32_t> cell_vertices = c_to_v->links(cell);
182 std::array<std::int32_t, 2> cols;
183 tcb::span<const std::int32_t> dofs1 = dofmap1->cell_dofs(cell);
184 for (std::int32_t i = 0; i < 2; ++i)
187 = std::find(cell_vertices.begin(), cell_vertices.end(), vertices[i]);
188 assert(it != cell_vertices.end());
189 const int local_vertex = std::distance(cell_vertices.begin(), it);
191 std::vector<std::int32_t>& local_v_dofs = local_vertex_dofs[local_vertex];
192 assert(local_v_dofs.size() == 1);
193 cols[i] = dofs1[local_v_dofs[0]];
196 if (global_indices[vertices[1]] < global_indices[vertices[0]])
201 mat_set(1, &row, 2, cols.data(), Ae.data());
This class represents a finite element function space defined by a mesh, a finite element,...
Definition: FunctionSpace.h:35
std::shared_ptr< const fem::DofMap > dofmap() const
The dofmap.
Definition: FunctionSpace.cpp:217
std::shared_ptr< const mesh::Mesh > mesh() const
The mesh.
Definition: FunctionSpace.cpp:210
This class provides a sparsity pattern data structure that can be used to initialize sparse matrices.
Definition: SparsityPattern.h:36
Finite element method functionality.
Definition: assemble_matrix_impl.h:23
void assemble_discrete_gradient(const std::function< int(std::int32_t, const std::int32_t *, std::int32_t, const std::int32_t *, const T *)> &mat_set, const fem::FunctionSpace &V0, const fem::FunctionSpace &V1)
Definition: discreteoperators.h:69
la::SparsityPattern create_sparsity_discrete_gradient(const fem::FunctionSpace &V0, const fem::FunctionSpace &V1)
Definition: discreteoperators.cpp:13
int cell_num_entities(mesh::CellType type, int dim)
Number of entities of dimension dim.
Definition: cell_types.cpp:194