DOLFINx
DOLFINx C++ interface
assemble_scalar_impl.h
1// Copyright (C) 2019-2020 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 "Constant.h"
10#include "Form.h"
11#include "FunctionSpace.h"
12#include "utils.h"
13#include <dolfinx/common/IndexMap.h>
14#include <dolfinx/common/utils.h>
15#include <dolfinx/mesh/Geometry.h>
16#include <dolfinx/mesh/Mesh.h>
17#include <dolfinx/mesh/Topology.h>
18#include <memory>
19#include <vector>
20
21namespace dolfinx::fem::impl
22{
23
25template <typename T>
26T assemble_cells(const mesh::Geometry& geometry,
27 const std::span<const std::int32_t>& cells,
28 const std::function<void(T*, const T*, const T*,
29 const scalar_value_type_t<T>*,
30 const int*, const std::uint8_t*)>& fn,
31 const std::span<const T>& constants,
32 const std::span<const T>& coeffs, int cstride)
33{
34 T value(0);
35 if (cells.empty())
36 return value;
37
38 // Prepare cell geometry
39 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
40 const std::size_t num_dofs_g = geometry.cmap().dim();
41 std::span<const double> x_g = geometry.x();
42
43 // Create data structures used in assembly
44 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
45
46 // Iterate over all cells
47 for (std::size_t index = 0; index < cells.size(); ++index)
48 {
49 std::int32_t c = cells[index];
50
51 // Get cell coordinates/geometry
52 auto x_dofs = x_dofmap.links(c);
53 for (std::size_t i = 0; i < x_dofs.size(); ++i)
54 {
55 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
56 std::next(coordinate_dofs.begin(), 3 * i));
57 }
58
59 const T* coeff_cell = coeffs.data() + index * cstride;
60 fn(&value, coeff_cell, constants.data(), coordinate_dofs.data(), nullptr,
61 nullptr);
62 }
63
64 return value;
65}
66
68template <typename T>
69T assemble_exterior_facets(
70 const mesh::Mesh& mesh, const std::span<const std::int32_t>& facets,
71 const std::function<void(T*, const T*, const T*,
72 const scalar_value_type_t<T>*, const int*,
73 const std::uint8_t*)>& fn,
74 const std::span<const T>& constants, const std::span<const T>& coeffs,
75 int cstride)
76{
77 T value(0);
78 if (facets.empty())
79 return value;
80
81 // Prepare cell geometry
82 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
83 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
84 std::span<const double> x_g = mesh.geometry().x();
85
86 // Create data structures used in assembly
87 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * num_dofs_g);
88
89 // Iterate over all facets
90 assert(facets.size() % 2 == 0);
91 for (std::size_t index = 0; index < facets.size(); index += 2)
92 {
93 std::int32_t cell = facets[index];
94 std::int32_t local_facet = facets[index + 1];
95
96 // Get cell coordinates/geometry
97 auto x_dofs = x_dofmap.links(cell);
98 for (std::size_t i = 0; i < x_dofs.size(); ++i)
99 {
100 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
101 std::next(coordinate_dofs.begin(), 3 * i));
102 }
103
104 const T* coeff_cell = coeffs.data() + index / 2 * cstride;
105 fn(&value, coeff_cell, constants.data(), coordinate_dofs.data(),
106 &local_facet, nullptr);
107 }
108
109 return value;
110}
111
113template <typename T>
114T assemble_interior_facets(
115 const mesh::Mesh& mesh, const std::span<const std::int32_t>& facets,
116 const std::function<void(T*, const T*, const T*,
117 const scalar_value_type_t<T>*, const int*,
118 const std::uint8_t*)>& fn,
119 const std::span<const T>& constants, const std::span<const T>& coeffs,
120 int cstride, const std::span<const int>& offsets,
121 const std::span<const std::uint8_t>& perms)
122{
123 T value(0);
124 if (facets.empty())
125 return value;
126
127 const int tdim = mesh.topology().dim();
128
129 // Prepare cell geometry
130 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
131 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
132 std::span<const double> x_g = mesh.geometry().x();
133
134 // Create data structures used in assembly
135 using X = scalar_value_type_t<T>;
136 std::vector<X> coordinate_dofs(2 * num_dofs_g * 3);
137 std::span<X> cdofs0(coordinate_dofs.data(), num_dofs_g * 3);
138 std::span<X> cdofs1(coordinate_dofs.data() + num_dofs_g * 3, num_dofs_g * 3);
139
140 std::vector<T> coeff_array(2 * offsets.back());
141 assert(offsets.back() == cstride);
142
143 const int num_cell_facets
144 = mesh::cell_num_entities(mesh.topology().cell_type(), tdim - 1);
145
146 // Iterate over all facets
147 assert(facets.size() % 4 == 0);
148 for (std::size_t index = 0; index < facets.size(); index += 4)
149 {
150 std::array<std::int32_t, 2> cells = {facets[index], facets[index + 2]};
151 std::array<std::int32_t, 2> local_facet
152 = {facets[index + 1], facets[index + 3]};
153
154 // Get cell geometry
155 auto x_dofs0 = x_dofmap.links(cells[0]);
156 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
157 {
158 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs0[i]),
159 std::next(cdofs0.begin(), 3 * i));
160 }
161 auto x_dofs1 = x_dofmap.links(cells[1]);
162 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
163 {
164 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs1[i]),
165 std::next(cdofs1.begin(), 3 * i));
166 }
167
168 const std::array perm{perms[cells[0] * num_cell_facets + local_facet[0]],
169 perms[cells[1] * num_cell_facets + local_facet[1]]};
170 fn(&value, coeffs.data() + index / 2 * cstride, constants.data(),
171 coordinate_dofs.data(), local_facet.data(), perm.data());
172 }
173
174 return value;
175}
176
178template <typename T>
180 const fem::Form<T>& M, const std::span<const T>& constants,
181 const std::map<std::pair<IntegralType, int>,
182 std::pair<std::span<const T>, int>>& coefficients)
183{
184 std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
185 assert(mesh);
186
187 T value = 0;
188 for (int i : M.integral_ids(IntegralType::cell))
189 {
190 const auto& fn = M.kernel(IntegralType::cell, i);
191 const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
192 const std::vector<std::int32_t>& cells = M.cell_domains(i);
193 value += impl::assemble_cells(mesh->geometry(), cells, fn, constants,
194 coeffs, cstride);
195 }
196
197 for (int i : M.integral_ids(IntegralType::exterior_facet))
198 {
199 const auto& fn = M.kernel(IntegralType::exterior_facet, i);
200 const auto& [coeffs, cstride]
201 = coefficients.at({IntegralType::exterior_facet, i});
202 const std::vector<std::int32_t>& facets = M.exterior_facet_domains(i);
203 value += impl::assemble_exterior_facets(*mesh, facets, fn, constants,
204 coeffs, cstride);
205 }
206
207 if (M.num_integrals(IntegralType::interior_facet) > 0)
208 {
209 mesh->topology_mutable().create_entity_permutations();
210
211 const std::vector<std::uint8_t>& perms
212 = mesh->topology().get_facet_permutations();
213
214 const std::vector<int> c_offsets = M.coefficient_offsets();
215 for (int i : M.integral_ids(IntegralType::interior_facet))
216 {
217 const auto& fn = M.kernel(IntegralType::interior_facet, i);
218 const auto& [coeffs, cstride]
219 = coefficients.at({IntegralType::interior_facet, i});
220 const std::vector<std::int32_t>& facets = M.interior_facet_domains(i);
221 value += impl::assemble_interior_facets(
222 *mesh, facets, fn, constants, coeffs, cstride, c_offsets, perms);
223 }
224 }
225
226 return value;
227}
228
229} // namespace dolfinx::fem::impl
void cells(la::SparsityPattern &pattern, const mesh::Topology &topology, const std::array< const std::reference_wrapper< const DofMap >, 2 > &dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:18
T assemble_scalar(const Form< T > &M, const std::span< const T > &constants, const std::map< std::pair< IntegralType, int >, std::pair< std::span< const T >, int > > &coefficients)
Assemble functional into scalar. The caller supplies the form constants and coefficients for this ver...
Definition: assembler.h:55
@ interior_facet
Interior facet.
@ exterior_facet
Exterior facet.
int cell_num_entities(CellType type, int dim)
Number of entities of dimension dim.
Definition: cell_types.cpp:182