Reference documentation for deal.II version 9.6.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
23
24
25namespace Particles
26{
27 namespace Utilities
28 {
29 template <int dim, int spacedim, typename number>
30 void
32 const DoFHandler<dim, spacedim> &space_dh,
33 const Particles::ParticleHandler<dim, spacedim> &particle_handler,
34 SparsityPatternBase &sparsity,
35 const AffineConstraints<number> &constraints,
36 const ComponentMask &space_comps)
37 {
38 if (particle_handler.n_locally_owned_particles() == 0)
39 return; // nothing to do here
40
41 const auto &fe = space_dh.get_fe();
42 const auto max_particles_per_cell =
43 particle_handler.n_global_max_particles_per_cell();
44
45 // Take care of components
46 const ComponentMask comps =
47 (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
48 space_comps);
49 AssertDimension(comps.size(), fe.n_components());
50
51 const auto n_comps = comps.n_selected_components();
52
53 // Global to local indices
54 std::vector<unsigned int> space_gtl(fe.n_components(),
56 for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
57 if (comps[i])
58 space_gtl[i] = j++;
59
60 // [TODO]: when the add_entries_local_to_global below will implement
61 // the version with the dof_mask, this should be uncommented.
62 // // Construct a dof_mask, used to distribute entries to the sparsity
63 // Table<2, bool> dof_mask(max_particles_per_cell * n_comps,
64 // fe.n_dofs_per_cell());
65 // dof_mask.fill(false);
66 // for (unsigned int i = 0; i < space_fe.n_dofs_per_cell(); ++i)
67 // {
68 // const auto comp_i = space_fe.system_to_component_index(i).first;
69 // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
70 // for (unsigned int j = 0; j < max_particles_per_cell; ++j)
71 // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true;
72 // }
73
74 std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
75 std::vector<types::particle_index> particle_indices(
76 max_particles_per_cell * n_comps);
77
78 auto particle = particle_handler.begin();
79 while (particle != particle_handler.end())
80 {
81 const auto &cell = particle->get_surrounding_cell();
82 const auto dh_cell =
83 typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
84 dh_cell->get_dof_indices(dof_indices);
85 const auto pic = particle_handler.particles_in_cell(cell);
86 const auto n_particles = particle_handler.n_particles_in_cell(cell);
87 particle_indices.resize(n_particles * n_comps);
88 Assert(pic.begin() == particle, ExcInternalError());
89 for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
90 {
91 const auto p_id = particle->get_id();
92 for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
93 {
94 const auto comp_j =
95 space_gtl[fe.system_to_component_index(j).first];
98 {p_id * n_comps + comp_j}, {dof_indices[j]}, sparsity);
99 }
100 }
101 // [TODO]: when this works, use this:
102 // constraints.add_entries_local_to_global(particle_indices,
103 // dof_indices,
104 // sparsity,
105 // dof_mask);
106 }
107 }
108
109
110
111 template <int dim, int spacedim, typename MatrixType>
112 void
114 const DoFHandler<dim, spacedim> &space_dh,
115 const Particles::ParticleHandler<dim, spacedim> &particle_handler,
116 MatrixType &matrix,
118 const ComponentMask &space_comps)
119 {
120 if (particle_handler.n_locally_owned_particles() == 0)
121 {
122 matrix.compress(VectorOperation::add);
123 return; // nothing else to do here
124 }
125
126 AssertDimension(matrix.n(), space_dh.n_dofs());
127
128 const auto &fe = space_dh.get_fe();
129 const auto max_particles_per_cell =
130 particle_handler.n_global_max_particles_per_cell();
131
132 // Take care of components
133 const ComponentMask comps =
134 (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
135 space_comps);
136 AssertDimension(comps.size(), fe.n_components());
137 const auto n_comps = comps.n_selected_components();
138
139 AssertDimension(matrix.m(),
140 particle_handler.n_global_particles() * n_comps);
141
142
143 // Global to local indices
144 std::vector<unsigned int> space_gtl(fe.n_components(),
146 for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
147 if (comps[i])
148 space_gtl[i] = j++;
149
150 // [TODO]: when the add_entries_local_to_global below will implement
151 // the version with the dof_mask, this should be uncommented.
152 // // Construct a dof_mask, used to distribute entries to the sparsity
153 // Table<2, bool> dof_mask(max_particles_per_cell * n_comps,
154 // fe.n_dofs_per_cell());
155 // dof_mask.fill(false);
156 // for (unsigned int i = 0; i < space_fe.n_dofs_per_cell(); ++i)
157 // {
158 // const auto comp_i = space_fe.system_to_component_index(i).first;
159 // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
160 // for (unsigned int j = 0; j < max_particles_per_cell; ++j)
161 // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true;
162 // }
163
164 std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
165 std::vector<types::particle_index> particle_indices(
166 max_particles_per_cell * n_comps);
167
169 max_particles_per_cell * n_comps, fe.n_dofs_per_cell());
170
171 auto particle = particle_handler.begin();
172 while (particle != particle_handler.end())
173 {
174 const auto &cell = particle->get_surrounding_cell();
175 const auto &dh_cell =
176 typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
177 dh_cell->get_dof_indices(dof_indices);
178 const auto pic = particle_handler.particles_in_cell(cell);
179 const auto n_particles = particle_handler.n_particles_in_cell(cell);
180 particle_indices.resize(n_particles * n_comps);
181 local_matrix.reinit({n_particles * n_comps, fe.n_dofs_per_cell()});
182 Assert(pic.begin() == particle, ExcInternalError());
183 for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
184 {
185 const auto &reference_location =
186 particle->get_reference_location();
187
188 for (unsigned int d = 0; d < n_comps; ++d)
189 particle_indices[i * n_comps + d] =
190 particle->get_id() * n_comps + d;
191
192 for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
193 {
194 const auto comp_j =
195 space_gtl[fe.system_to_component_index(j).first];
196 if (comp_j != numbers::invalid_unsigned_int)
197 local_matrix(i * n_comps + comp_j, j) =
198 fe.shape_value(j, reference_location);
199 }
200 }
201 constraints.distribute_local_to_global(local_matrix,
202 particle_indices,
203 dof_indices,
204 matrix);
205 }
206 matrix.compress(VectorOperation::add);
207 }
208
209#include "utilities.inst"
210
211 } // namespace Utilities
212} // 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:503
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:504
#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:113
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:31
static const unsigned int invalid_unsigned_int
Definition types.h:220