Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3110-g10dd77059b 2025-04-22 10:30:00+00:00
\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}} \newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=} \newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]} \newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
utilities.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#include <deal.II/base/config.h>
16
18
20
21#include <vector>
22
23
25
26
27namespace Particles
28{
29 namespace Utilities
30 {
31 template <int dim, int spacedim, typename number>
32 void
34 const DoFHandler<dim, spacedim> &space_dh,
35 const Particles::ParticleHandler<dim, spacedim> &particle_handler,
36 SparsityPatternBase &sparsity,
37 const AffineConstraints<number> &constraints,
38 const ComponentMask &space_comps)
39 {
40 if (particle_handler.n_locally_owned_particles() == 0)
41 return; // nothing to do here
42
43 const auto &fe = space_dh.get_fe();
44 const auto max_particles_per_cell =
45 particle_handler.n_global_max_particles_per_cell();
46
47 // Take care of components
48 const ComponentMask comps =
49 (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
50 space_comps);
51 AssertDimension(comps.size(), fe.n_components());
52
53 const auto n_comps = comps.n_selected_components();
54
55 // Global to local indices
56 std::vector<unsigned int> space_gtl(fe.n_components(),
58 for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
59 if (comps[i])
60 space_gtl[i] = j++;
61
62 // [TODO]: when the add_entries_local_to_global below will implement
63 // the version with the dof_mask, this should be uncommented.
64 // // Construct a dof_mask, used to distribute entries to the sparsity
65 // Table<2, bool> dof_mask(max_particles_per_cell * n_comps,
66 // fe.n_dofs_per_cell());
67 // dof_mask.fill(false);
68 // for (unsigned int i = 0; i < space_fe.n_dofs_per_cell(); ++i)
69 // {
70 // const auto comp_i = space_fe.system_to_component_index(i).first;
71 // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
72 // for (unsigned int j = 0; j < max_particles_per_cell; ++j)
73 // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true;
74 // }
75
76 std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
77 std::vector<types::particle_index> particle_indices(
78 max_particles_per_cell * n_comps);
79
80 auto particle = particle_handler.begin();
81 while (particle != particle_handler.end())
82 {
83 const auto &cell = particle->get_surrounding_cell();
84 const auto dh_cell =
85 typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
86 dh_cell->get_dof_indices(dof_indices);
87 const auto pic = particle_handler.particles_in_cell(cell);
88 const auto n_particles = particle_handler.n_particles_in_cell(cell);
89 particle_indices.resize(n_particles * n_comps);
90 Assert(pic.begin() == particle, ExcInternalError());
91 for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
92 {
93 const auto p_id = particle->get_id();
94 for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
95 {
96 const auto comp_j =
97 space_gtl[fe.system_to_component_index(j).first];
100 {p_id * n_comps + comp_j}, {dof_indices[j]}, sparsity);
101 }
102 }
103 // [TODO]: when this works, use this:
104 // constraints.add_entries_local_to_global(particle_indices,
105 // dof_indices,
106 // sparsity,
107 // dof_mask);
108 }
109 }
110
111
112
113 template <int dim, int spacedim, typename MatrixType>
114 void
116 const DoFHandler<dim, spacedim> &space_dh,
117 const Particles::ParticleHandler<dim, spacedim> &particle_handler,
118 MatrixType &matrix,
120 const ComponentMask &space_comps)
121 {
122 if (particle_handler.n_locally_owned_particles() == 0)
123 {
124 matrix.compress(VectorOperation::add);
125 return; // nothing else to do here
126 }
127
128 AssertDimension(matrix.n(), space_dh.n_dofs());
129
130 const auto &fe = space_dh.get_fe();
131 const auto max_particles_per_cell =
132 particle_handler.n_global_max_particles_per_cell();
133
134 // Take care of components
135 const ComponentMask comps =
136 (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
137 space_comps);
138 AssertDimension(comps.size(), fe.n_components());
139 const auto n_comps = comps.n_selected_components();
140
141 AssertDimension(matrix.m(),
142 particle_handler.n_global_particles() * n_comps);
143
144
145 // Global to local indices
146 std::vector<unsigned int> space_gtl(fe.n_components(),
148 for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
149 if (comps[i])
150 space_gtl[i] = j++;
151
152 // [TODO]: when the add_entries_local_to_global below will implement
153 // the version with the dof_mask, this should be uncommented.
154 // // Construct a dof_mask, used to distribute entries to the sparsity
155 // Table<2, bool> dof_mask(max_particles_per_cell * n_comps,
156 // fe.n_dofs_per_cell());
157 // dof_mask.fill(false);
158 // for (unsigned int i = 0; i < space_fe.n_dofs_per_cell(); ++i)
159 // {
160 // const auto comp_i = space_fe.system_to_component_index(i).first;
161 // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
162 // for (unsigned int j = 0; j < max_particles_per_cell; ++j)
163 // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true;
164 // }
165
166 std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
167 std::vector<types::particle_index> particle_indices(
168 max_particles_per_cell * n_comps);
169
171 max_particles_per_cell * n_comps, fe.n_dofs_per_cell());
172
173 auto particle = particle_handler.begin();
174 while (particle != particle_handler.end())
175 {
176 const auto &cell = particle->get_surrounding_cell();
177 const auto &dh_cell =
178 typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
179 dh_cell->get_dof_indices(dof_indices);
180 const auto pic = particle_handler.particles_in_cell(cell);
181 const auto n_particles = particle_handler.n_particles_in_cell(cell);
182 particle_indices.resize(n_particles * n_comps);
183 local_matrix.reinit({n_particles * n_comps, fe.n_dofs_per_cell()});
184 Assert(pic.begin() == particle, ExcInternalError());
185 for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
186 {
187 const auto &reference_location =
188 particle->get_reference_location();
189
190 for (unsigned int d = 0; d < n_comps; ++d)
191 particle_indices[i * n_comps + d] =
192 particle->get_id() * n_comps + d;
193
194 for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
195 {
196 const auto comp_j =
197 space_gtl[fe.system_to_component_index(j).first];
198 if (comp_j != numbers::invalid_unsigned_int)
199 local_matrix(i * n_comps + comp_j, j) =
200 fe.shape_value(j, reference_location);
201 }
202 }
203 constraints.distribute_local_to_global(local_matrix,
204 particle_indices,
205 dof_indices,
206 matrix);
207 }
208 matrix.compress(VectorOperation::add);
209 }
210
211#include "particles/utilities.inst"
212
213 } // namespace Utilities
214} // namespace Particles
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternBase &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
unsigned int size() const
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
const FiniteElement< dim, spacedim > & get_fe(const types::fe_index index=0) const
types::global_dof_index n_dofs() const
types::particle_index n_global_particles() const
particle_iterator begin() const
particle_iterator end() const
types::particle_index n_locally_owned_particles() const
types::particle_index n_particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell) const
particle_iterator_range particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
types::particle_index n_global_max_particles_per_cell() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
static ::ExceptionBase & ExcInternalError()
typename ActiveSelector::cell_iterator cell_iterator
void create_interpolation_matrix(const DoFHandler< dim, spacedim > &space_dh, const Particles::ParticleHandler< dim, spacedim > &particle_handler, MatrixType &matrix, const AffineConstraints< typename MatrixType::value_type > &constraints=AffineConstraints< typename MatrixType::value_type >(), const ComponentMask &space_comps={})
Definition utilities.cc:115
void create_interpolation_sparsity_pattern(const DoFHandler< dim, spacedim > &space_dh, const Particles::ParticleHandler< dim, spacedim > &particle_handler, SparsityPatternBase &sparsity, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const ComponentMask &space_comps={})
Definition utilities.cc:33
constexpr unsigned int invalid_unsigned_int
Definition types.h:238