DOLFINx 0.9.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
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 <algorithm>
14#include <dolfinx/common/IndexMap.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 <dolfinx::scalar T>
26T assemble_cells(mdspan2_t x_dofmap, std::span<const scalar_value_type_t<T>> x,
27 std::span<const std::int32_t> cells, FEkernel<T> auto fn,
28 std::span<const T> constants, std::span<const T> coeffs,
29 int cstride)
30{
31 T value(0);
32 if (cells.empty())
33 return value;
34
35 // Create data structures used in assembly
36 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * x_dofmap.extent(1));
37
38 // Iterate over all cells
39 for (std::size_t index = 0; index < cells.size(); ++index)
40 {
41 std::int32_t c = cells[index];
42
43 // Get cell coordinates/geometry
44 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
45 x_dofmap, c, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
46 for (std::size_t i = 0; i < x_dofs.size(); ++i)
47 {
48 std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
49 std::next(coordinate_dofs.begin(), 3 * i));
50 }
51
52 const T* coeff_cell = coeffs.data() + index * cstride;
53 fn(&value, coeff_cell, constants.data(), coordinate_dofs.data(), nullptr,
54 nullptr);
55 }
56
57 return value;
58}
59
61template <dolfinx::scalar T>
62T assemble_exterior_facets(mdspan2_t x_dofmap,
63 std::span<const scalar_value_type_t<T>> x,
64 int num_facets_per_cell,
65 std::span<const std::int32_t> facets,
66 FEkernel<T> auto fn, std::span<const T> constants,
67 std::span<const T> coeffs, int cstride,
68 std::span<const std::uint8_t> perms)
69{
70 T value(0);
71 if (facets.empty())
72 return value;
73
74 // Create data structures used in assembly
75 std::vector<scalar_value_type_t<T>> coordinate_dofs(3 * x_dofmap.extent(1));
76
77 // Iterate over all facets
78 assert(facets.size() % 2 == 0);
79 for (std::size_t index = 0; index < facets.size(); index += 2)
80 {
81 std::int32_t cell = facets[index];
82 std::int32_t local_facet = facets[index + 1];
83
84 // Get cell coordinates/geometry
85 auto x_dofs = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
86 x_dofmap, cell, MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
87 for (std::size_t i = 0; i < x_dofs.size(); ++i)
88 {
89 std::copy_n(std::next(x.begin(), 3 * x_dofs[i]), 3,
90 std::next(coordinate_dofs.begin(), 3 * i));
91 }
92
93 // Permutations
94 std::uint8_t perm
95 = perms.empty() ? 0 : perms[cell * num_facets_per_cell + local_facet];
96 const T* coeff_cell = coeffs.data() + index / 2 * cstride;
97 fn(&value, coeff_cell, constants.data(), coordinate_dofs.data(),
98 &local_facet, &perm);
99 }
100
101 return value;
102}
103
105template <dolfinx::scalar T>
106T assemble_interior_facets(mdspan2_t x_dofmap,
107 std::span<const scalar_value_type_t<T>> x,
108 int num_facets_per_cell,
109 std::span<const std::int32_t> facets,
110 FEkernel<T> auto fn, std::span<const T> constants,
111 std::span<const T> coeffs, int cstride,
112 std::span<const int> offsets,
113 std::span<const std::uint8_t> perms)
114{
115 T value(0);
116 if (facets.empty())
117 return value;
118
119 // Create data structures used in assembly
120 using X = scalar_value_type_t<T>;
121 std::vector<X> coordinate_dofs(2 * x_dofmap.extent(1) * 3);
122 std::span<X> cdofs0(coordinate_dofs.data(), x_dofmap.extent(1) * 3);
123 std::span<X> cdofs1(coordinate_dofs.data() + x_dofmap.extent(1) * 3,
124 x_dofmap.extent(1) * 3);
125
126 std::vector<T> coeff_array(2 * offsets.back());
127 assert(offsets.back() == cstride);
128
129 // Iterate over all facets
130 assert(facets.size() % 4 == 0);
131 for (std::size_t index = 0; index < facets.size(); index += 4)
132 {
133 std::array<std::int32_t, 2> cells = {facets[index], facets[index + 2]};
134 std::array<std::int32_t, 2> local_facet
135 = {facets[index + 1], facets[index + 3]};
136
137 // Get cell geometry
138 auto x_dofs0 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
139 x_dofmap, cells[0], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
140 for (std::size_t i = 0; i < x_dofs0.size(); ++i)
141 {
142 std::copy_n(std::next(x.begin(), 3 * x_dofs0[i]), 3,
143 std::next(cdofs0.begin(), 3 * i));
144 }
145 auto x_dofs1 = MDSPAN_IMPL_STANDARD_NAMESPACE::submdspan(
146 x_dofmap, cells[1], MDSPAN_IMPL_STANDARD_NAMESPACE::full_extent);
147 for (std::size_t i = 0; i < x_dofs1.size(); ++i)
148 {
149 std::copy_n(std::next(x.begin(), 3 * x_dofs1[i]), 3,
150 std::next(cdofs1.begin(), 3 * i));
151 }
152
153 std::array perm
154 = perms.empty()
155 ? std::array<std::uint8_t, 2>{0, 0}
156 : std::array{
157 perms[cells[0] * num_facets_per_cell + local_facet[0]],
158 perms[cells[1] * num_facets_per_cell + local_facet[1]]};
159 fn(&value, coeffs.data() + index / 2 * cstride, constants.data(),
160 coordinate_dofs.data(), local_facet.data(), perm.data());
161 }
162
163 return value;
164}
165
167template <dolfinx::scalar T, std::floating_point U>
168T assemble_scalar(
169 const fem::Form<T, U>& M, mdspan2_t x_dofmap,
170 std::span<const scalar_value_type_t<T>> x, std::span<const T> constants,
171 const std::map<std::pair<IntegralType, int>,
172 std::pair<std::span<const T>, int>>& coefficients)
173{
174 std::shared_ptr<const mesh::Mesh<U>> mesh = M.mesh();
175 assert(mesh);
176
177 T value = 0;
178 for (int i : M.integral_ids(IntegralType::cell))
179 {
180 auto fn = M.kernel(IntegralType::cell, i);
181 assert(fn);
182 auto& [coeffs, cstride] = coefficients.at({IntegralType::cell, i});
183 std::span<const std::int32_t> cells = M.domain(IntegralType::cell, i);
184 value += impl::assemble_cells(x_dofmap, x, cells, fn, constants, coeffs,
185 cstride);
186 }
187
188 std::span<const std::uint8_t> perms;
189 if (M.needs_facet_permutations())
190 {
191 mesh->topology_mutable()->create_entity_permutations();
192 perms = std::span(mesh->topology()->get_facet_permutations());
193 }
194
195 mesh::CellType cell_type = mesh->topology()->cell_type();
196 int num_facets_per_cell
197 = mesh::cell_num_entities(cell_type, mesh->topology()->dim() - 1);
198 for (int i : M.integral_ids(IntegralType::exterior_facet))
199 {
200 auto fn = M.kernel(IntegralType::exterior_facet, i);
201 assert(fn);
202 auto& [coeffs, cstride]
203 = coefficients.at({IntegralType::exterior_facet, i});
204 value += impl::assemble_exterior_facets(
205 x_dofmap, x, num_facets_per_cell,
206 M.domain(IntegralType::exterior_facet, i), fn, constants, coeffs,
207 cstride, perms);
208 }
209
210 for (int i : M.integral_ids(IntegralType::interior_facet))
211 {
212 const std::vector<int> c_offsets = M.coefficient_offsets();
213 auto fn = M.kernel(IntegralType::interior_facet, i);
214 assert(fn);
215 auto& [coeffs, cstride]
216 = coefficients.at({IntegralType::interior_facet, i});
217 value += impl::assemble_interior_facets(
218 x_dofmap, x, num_facets_per_cell,
219 M.domain(IntegralType::interior_facet, i), fn, constants, coeffs,
220 cstride, c_offsets, perms);
221 }
222
223 return value;
224}
225
226} // namespace dolfinx::fem::impl
Functions supporting finite element method operations.
void cells(la::SparsityPattern &pattern, std::array< std::span< const std::int32_t >, 2 > cells, std::array< std::reference_wrapper< const DofMap >, 2 > dofmaps)
Iterate over cells and insert entries into sparsity pattern.
Definition sparsitybuild.cpp:16
@ interior_facet
Interior facet.
@ exterior_facet
Exterior facet.
int cell_num_entities(CellType type, int dim)
Definition cell_types.cpp:139
CellType
Cell type identifier.
Definition cell_types.h:22