DOLFIN-X
DOLFIN-X C++ interface
discreteoperators.h
1 // Copyright (C) 2015 Garth N. Wells
2 //
3 // This file is part of DOLFINX (https://www.fenicsproject.org)
4 //
5 // SPDX-License-Identifier: LGPL-3.0-or-later
6 
7 #pragma once
8 
9 #include <array>
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>
15 #include <memory>
16 #include <vector>
17 
18 namespace dolfinx::fem
19 {
20 
39  const fem::FunctionSpace& V1);
40 
60 template <typename T>
62  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
63  const std::int32_t*, const T*)>& mat_set,
64  const fem::FunctionSpace& V0, const fem::FunctionSpace& V1);
65 } // namespace dolfinx::fem
66 
67 using namespace dolfinx;
68 template <typename T>
70  const std::function<int(std::int32_t, const std::int32_t*, std::int32_t,
71  const std::int32_t*, const T*)>& mat_set,
72  const fem::FunctionSpace& V0, const fem::FunctionSpace& V1)
73 {
74  // Get mesh
75  std::shared_ptr<const mesh::Mesh> mesh = V0.mesh();
76  assert(mesh);
77 
78  // Check that mesh is the same for both function spaces
79  assert(V1.mesh());
80  if (mesh != V1.mesh())
81  {
82  throw std::runtime_error("Compute discrete gradient operator. Function "
83  "spaces do not share the same mesh");
84  }
85 
86  // Check that V0 is a (lowest-order) edge basis
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)
92  {
93  throw std::runtime_error(
94  "Cannot compute discrete gradient operator. Function "
95  "spaces is not a lowest-order edge space");
96  }
97 
98  // Check that V1 is a linear nodal basis
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)
104  {
105  throw std::runtime_error(
106  "Cannot compute discrete gradient operator. Function "
107  "space is not a linear nodal function space");
108  }
109 
110  // Build maps from entities to local dof indices
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;
115 
116  // Copy index maps from dofmaps
117  std::array<std::shared_ptr<const common::IndexMap>, 2> index_maps
118  = {{V0.dofmap()->index_map, V1.dofmap()->index_map}};
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]);
124 
125  // Initialize required connectivities
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);
135 
136  // Build sparsity pattern
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();
140  assert(dofmap0);
141  // Create local lookup table for local edge to cell dofs
142  const int num_edges_per_cell
143  = mesh::cell_num_entities(mesh->topology().cell_type(), 1);
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);
147  // Create local lookup table for local vertex to cell dofs
148  const int num_vertices_per_cell
149  = mesh::cell_num_entities(mesh->topology().cell_type(), 0);
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);
153 
154  // Build discrete gradient operator/matrix
155  const std::shared_ptr<const fem::DofMap> dofmap1 = V1.dofmap();
156  assert(dofmap1);
157  const std::vector<std::int64_t>& global_indices
158  = mesh->topology().index_map(0)->global_indices();
159  std::array<T, 2> Ae;
160  for (std::int32_t e = 0; e < num_edges; ++e)
161  {
162  // Find local index of edge in one of the cells it is part of
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);
170 
171  // Find the dofs located on the edge
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]];
176 
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);
180 
181  // Find local index of each of the vertices and map to local dof
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)
185  {
186  const auto it
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);
190 
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]];
194  }
195 
196  if (global_indices[vertices[1]] < global_indices[vertices[0]])
197  Ae = {1, -1};
198  else
199  Ae = {-1, 1};
200 
201  mat_set(1, &row, 2, cols.data(), Ae.data());
202  }
203 }
204 //-----------------------------------------------------------------------------
This class represents a finite element function space defined by a mesh, a finite element,...
Definition: FunctionSpace.h:34
std::shared_ptr< const fem::DofMap > dofmap() const
The dofmap.
Definition: FunctionSpace.cpp:202
std::shared_ptr< const mesh::Mesh > mesh() const
The mesh.
Definition: FunctionSpace.cpp:195
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:182