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 xtl::span<const std::int32_t>& cells,
28 const std::function<void(T*, const T*, const T*, const double*,
29 const int*, const std::uint8_t*)>& fn,
30 const xtl::span<const T>& constants,
31 const xtl::span<const T>& coeffs, int cstride)
32{
33 T value(0);
34 if (cells.empty())
35 return value;
36
37 // Prepare cell geometry
38 const graph::AdjacencyList<std::int32_t>& x_dofmap = geometry.dofmap();
39 const std::size_t num_dofs_g = geometry.cmap().dim();
40 xtl::span<const double> x_g = geometry.x();
41
42 // Create data structures used in assembly
43 std::vector<double> coordinate_dofs(3 * num_dofs_g);
44
45 // Iterate over all cells
46 for (std::size_t index = 0; index < cells.size(); ++index)
47 {
48 std::int32_t c = cells[index];
49
50 // Get cell coordinates/geometry
51 auto x_dofs = x_dofmap.links(c);
52 for (std::size_t i = 0; i < x_dofs.size(); ++i)
53 {
54 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
55 std::next(coordinate_dofs.begin(), 3 * i));
56 }
57
58 const T* coeff_cell = coeffs.data() + index * cstride;
59 fn(&value, coeff_cell, constants.data(), coordinate_dofs.data(), nullptr,
60 nullptr);
61 }
62
63 return value;
64}
65
67template <typename T>
68T assemble_exterior_facets(
69 const mesh::Mesh& mesh,
70 const xtl::span<const std::pair<std::int32_t, int>>& facets,
71 const std::function<void(T*, const T*, const T*, const double*, const int*,
72 const std::uint8_t*)>& fn,
73 const xtl::span<const T>& constants, const xtl::span<const T>& coeffs,
74 int cstride)
75{
76 T value(0);
77 if (facets.empty())
78 return value;
79
80 // Prepare cell geometry
81 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
82 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
83 xtl::span<const double> x_g = mesh.geometry().x();
84
85 // Create data structures used in assembly
86 std::vector<double> coordinate_dofs(3 * num_dofs_g);
87
88 // Iterate over all facets
89 for (std::size_t index = 0; index < facets.size(); ++index)
90 {
91 std::int32_t cell = facets[index].first;
92 int local_facet = facets[index].second;
93
94 // Get cell coordinates/geometry
95 auto x_dofs = x_dofmap.links(cell);
96 for (std::size_t i = 0; i < x_dofs.size(); ++i)
97 {
98 common::impl::copy_N<3>(std::next(x_g.begin(), 3 * x_dofs[i]),
99 std::next(coordinate_dofs.begin(), 3 * i));
100 }
101
102 const T* coeff_cell = coeffs.data() + index * cstride;
103 fn(&value, coeff_cell, constants.data(), coordinate_dofs.data(),
104 &local_facet, nullptr);
105 }
106
107 return value;
108}
109
111template <typename T>
112T assemble_interior_facets(
113 const mesh::Mesh& mesh,
114 const xtl::span<const std::tuple<std::int32_t, int, std::int32_t, int>>&
115 facets,
116 const std::function<void(T*, const T*, const T*, const double*, const int*,
117 const std::uint8_t*)>& fn,
118 const xtl::span<const T>& constants, const xtl::span<const T>& coeffs,
119 int cstride, const xtl::span<const int>& offsets,
120 const xtl::span<const std::uint8_t>& perms)
121{
122 T value(0);
123 if (facets.empty())
124 return value;
125
126 const int tdim = mesh.topology().dim();
127
128 // Prepare cell geometry
129 const graph::AdjacencyList<std::int32_t>& x_dofmap = mesh.geometry().dofmap();
130 const std::size_t num_dofs_g = mesh.geometry().cmap().dim();
131 xtl::span<const double> x_g = mesh.geometry().x();
132
133 // Create data structures used in assembly
134 xt::xtensor<double, 3> coordinate_dofs({2, num_dofs_g, 3});
135 std::vector<T> coeff_array(2 * offsets.back());
136 assert(offsets.back() == cstride);
137
138 const int num_cell_facets
139 = mesh::cell_num_entities(mesh.topology().cell_type(), tdim - 1);
140
141 // Iterate over all facets
142 for (std::size_t index = 0; index < facets.size(); ++index)
143 {
144 const std::array<std::int32_t, 2> cells
145 = {std::get<0>(facets[index]), std::get<2>(facets[index])};
146 const std::array<int, 2> local_facet
147 = {std::get<1>(facets[index]), std::get<3>(facets[index])};
148
149 // Get cell geometry
150 auto x_dofs0 = x_dofmap.links(cells[0]);
151 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
152 {
153 common::impl::copy_N<3>(
154 std::next(x_g.begin(), 3 * x_dofs0[i]),
155 xt::view(coordinate_dofs, 0, i, xt::all()).begin());
156 }
157 auto x_dofs1 = x_dofmap.links(cells[1]);
158 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
159 {
160 common::impl::copy_N<3>(
161 std::next(x_g.begin(), 3 * x_dofs1[i]),
162 xt::view(coordinate_dofs, 1, i, xt::all()).begin());
163 }
164
165 const std::array perm{perms[cells[0] * num_cell_facets + local_facet[0]],
166 perms[cells[1] * num_cell_facets + local_facet[1]]};
167 fn(&value, coeffs.data() + index * 2 * cstride, constants.data(),
168 coordinate_dofs.data(), local_facet.data(), perm.data());
169 }
170
171 return value;
172}
173
175template <typename T>
177 const fem::Form<T>& M, const xtl::span<const T>& constants,
178 const std::map<std::pair<IntegralType, int>,
179 std::pair<xtl::span<const T>, int>>& coefficients)
180{
181 std::shared_ptr<const mesh::Mesh> mesh = M.mesh();
182 assert(mesh);
183
184 T value = 0;
185 for (int i : M.integral_ids(IntegralType::cell))
186 {
187 const auto& fn = M.kernel(IntegralType::cell, i);
188 const auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
189 const std::vector<std::int32_t>& cells = M.cell_domains(i);
190 value += impl::assemble_cells(mesh->geometry(), cells, fn, constants,
191 coeffs, cstride);
192 }
193
194 for (int i : M.integral_ids(IntegralType::exterior_facet))
195 {
196 const auto& fn = M.kernel(IntegralType::exterior_facet, i);
197 const auto& [coeffs, cstride]
198 = coefficients.at({IntegralType::exterior_facet, i});
199 const std::vector<std::pair<std::int32_t, int>>& facets
200 = M.exterior_facet_domains(i);
201 value += impl::assemble_exterior_facets(*mesh, facets, fn, constants,
202 coeffs, cstride);
203 }
204
205 if (M.num_integrals(IntegralType::interior_facet) > 0)
206 {
207 mesh->topology_mutable().create_entity_permutations();
208
209 const std::vector<std::uint8_t>& perms
210 = mesh->topology().get_facet_permutations();
211
212 const std::vector<int> c_offsets = M.coefficient_offsets();
213 for (int i : M.integral_ids(IntegralType::interior_facet))
214 {
215 const auto& fn = M.kernel(IntegralType::interior_facet, i);
216 const auto& [coeffs, cstride]
217 = coefficients.at({IntegralType::interior_facet, i});
218 const std::vector<std::tuple<std::int32_t, int, std::int32_t, int>>&
219 facets
220 = 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 fem::DofMap >, 2 > &dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition: sparsitybuild.cpp:18
T assemble_scalar(const Form< T > &M, const xtl::span< const T > &constants, const std::map< std::pair< IntegralType, int >, std::pair< xtl::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:185