Reference documentation for deal.II version Git d745029578 2021-12-01 09:26:47 +0100
\(\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\}}\)
utilities.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2020 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/config.h>
17 
19 
21 
22 
24 
25 
26 namespace Particles
27 {
28  namespace Utilities
29  {
30  template <int dim, int spacedim, typename SparsityType, typename number>
31  void
33  const DoFHandler<dim, spacedim> & space_dh,
34  const Particles::ParticleHandler<dim, spacedim> &particle_handler,
35  SparsityType & sparsity,
36  const AffineConstraints<number> & constraints,
37  const ComponentMask & space_comps)
38  {
39  if (particle_handler.n_locally_owned_particles() == 0)
40  return; // nothing to do here
41 
42  const auto &fe = space_dh.get_fe();
43  const auto max_particles_per_cell =
44  particle_handler.n_global_max_particles_per_cell();
45 
46  // Take care of components
47  const ComponentMask comps =
48  (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
49  space_comps);
50  AssertDimension(comps.size(), fe.n_components());
51 
52  const auto n_comps = comps.n_selected_components();
53 
54  // Global to local indices
55  std::vector<unsigned int> space_gtl(fe.n_components(),
57  for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
58  if (comps[i])
59  space_gtl[i] = j++;
60 
61  // [TODO]: when the add_entries_local_to_global below will implement
62  // the version with the dof_mask, this should be uncommented.
63  // // Construct a dof_mask, used to distribute entries to the sparsity
64  // Table<2, bool> dof_mask(max_particles_per_cell * n_comps,
65  // fe.n_dofs_per_cell());
66  // dof_mask.fill(false);
67  // for (unsigned int i = 0; i < space_fe.n_dofs_per_cell(); ++i)
68  // {
69  // const auto comp_i = space_fe.system_to_component_index(i).first;
70  // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
71  // for (unsigned int j = 0; j < max_particles_per_cell; ++j)
72  // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true;
73  // }
74 
75  std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
76  std::vector<types::particle_index> particle_indices(
77  max_particles_per_cell * n_comps);
78 
79  auto particle = particle_handler.begin();
80  while (particle != particle_handler.end())
81  {
82  const auto &cell = particle->get_surrounding_cell();
83  const auto dh_cell =
84  typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
85  dh_cell->get_dof_indices(dof_indices);
86  const auto pic = particle_handler.particles_in_cell(cell);
87  const auto n_particles = particle_handler.n_particles_in_cell(cell);
88  particle_indices.resize(n_particles * n_comps);
89  Assert(pic.begin() == particle, ExcInternalError());
90  for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
91  {
92  const auto p_id = particle->get_id();
93  for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
94  {
95  const auto comp_j =
96  space_gtl[fe.system_to_component_index(j).first];
97  if (comp_j != numbers::invalid_unsigned_int)
98  constraints.add_entries_local_to_global(
99  {p_id * n_comps + comp_j}, {dof_indices[j]}, sparsity);
100  }
101  }
102  // [TODO]: when this works, use this:
103  // constraints.add_entries_local_to_global(particle_indices,
104  // dof_indices,
105  // sparsity,
106  // dof_mask);
107  }
108  }
109 
110 
111 
112  template <int dim, int spacedim, typename MatrixType>
113  void
115  const DoFHandler<dim, spacedim> & space_dh,
116  const Particles::ParticleHandler<dim, spacedim> &particle_handler,
117  MatrixType & matrix,
119  const ComponentMask & space_comps)
120  {
121  if (particle_handler.n_locally_owned_particles() == 0)
122  {
123  matrix.compress(VectorOperation::add);
124  return; // nothing else to do here
125  }
126 
127  AssertDimension(matrix.n(), space_dh.n_dofs());
128 
129  const auto &fe = space_dh.get_fe();
130  const auto max_particles_per_cell =
131  particle_handler.n_global_max_particles_per_cell();
132 
133  // Take care of components
134  const ComponentMask comps =
135  (space_comps.size() == 0 ? ComponentMask(fe.n_components(), true) :
136  space_comps);
137  AssertDimension(comps.size(), fe.n_components());
138  const auto n_comps = comps.n_selected_components();
139 
140  AssertDimension(matrix.m(),
141  particle_handler.n_global_particles() * n_comps);
142 
143 
144  // Global to local indices
145  std::vector<unsigned int> space_gtl(fe.n_components(),
147  for (unsigned int i = 0, j = 0; i < space_gtl.size(); ++i)
148  if (comps[i])
149  space_gtl[i] = j++;
150 
151  // [TODO]: when the add_entries_local_to_global below will implement
152  // the version with the dof_mask, this should be uncommented.
153  // // Construct a dof_mask, used to distribute entries to the sparsity
154  // Table<2, bool> dof_mask(max_particles_per_cell * n_comps,
155  // fe.n_dofs_per_cell());
156  // dof_mask.fill(false);
157  // for (unsigned int i = 0; i < space_fe.n_dofs_per_cell(); ++i)
158  // {
159  // const auto comp_i = space_fe.system_to_component_index(i).first;
160  // if (space_gtl[comp_i] != numbers::invalid_unsigned_int)
161  // for (unsigned int j = 0; j < max_particles_per_cell; ++j)
162  // dof_mask(i, j * n_comps + space_gtl[comp_i]) = true;
163  // }
164 
165  std::vector<types::global_dof_index> dof_indices(fe.n_dofs_per_cell());
166  std::vector<types::particle_index> particle_indices(
167  max_particles_per_cell * n_comps);
168 
170  max_particles_per_cell * n_comps, fe.n_dofs_per_cell());
171 
172  auto particle = particle_handler.begin();
173  while (particle != particle_handler.end())
174  {
175  const auto &cell = particle->get_surrounding_cell();
176  const auto &dh_cell =
177  typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
178  dh_cell->get_dof_indices(dof_indices);
179  const auto pic = particle_handler.particles_in_cell(cell);
180  const auto n_particles = particle_handler.n_particles_in_cell(cell);
181  particle_indices.resize(n_particles * n_comps);
182  local_matrix.reinit({n_particles * n_comps, fe.n_dofs_per_cell()});
183  Assert(pic.begin() == particle, ExcInternalError());
184  for (unsigned int i = 0; particle != pic.end(); ++particle, ++i)
185  {
186  const auto &reference_location =
187  particle->get_reference_location();
188 
189  for (unsigned int d = 0; d < n_comps; ++d)
190  particle_indices[i * n_comps + d] =
191  particle->get_id() * n_comps + d;
192 
193  for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j)
194  {
195  const auto comp_j =
196  space_gtl[fe.system_to_component_index(j).first];
197  if (comp_j != numbers::invalid_unsigned_int)
198  local_matrix(i * n_comps + comp_j, j) =
199  fe.shape_value(j, reference_location);
200  }
201  }
202  constraints.distribute_local_to_global(local_matrix,
203  particle_indices,
204  dof_indices,
205  matrix);
206  }
207  matrix.compress(VectorOperation::add);
208  }
209 
210 #include "utilities.inst"
211 
212  } // namespace Utilities
213 } // namespace Particles
types::particle_index n_global_particles() const
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
types::particle_index n_locally_owned_particles() const
Contents is actually a matrix.
types::particle_index n_global_max_particles_per_cell() const
types::particle_index n_particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell) const
void distribute_local_to_global(const InVector &local_vector, const std::vector< size_type > &local_dof_indices, OutVector &global_vector) const
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=ComponentMask())
Definition: utilities.cc:114
const FiniteElement< dim, spacedim > & get_fe(const unsigned int index=0) const
#define Assert(cond, exc)
Definition: exceptions.h:1461
particle_iterator_range particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:404
types::global_dof_index n_dofs() const
particle_iterator begin() const
unsigned int size() const
iterator end(const size_type r)
void add_entries_local_to_global(const std::vector< size_type > &local_dof_indices, SparsityPatternType &sparsity_pattern, const bool keep_constrained_entries=true, const Table< 2, bool > &dof_mask=Table< 2, bool >()) const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:403
unsigned int n_selected_components(const unsigned int overall_number_of_components=numbers::invalid_unsigned_int) const
void create_interpolation_sparsity_pattern(const DoFHandler< dim, spacedim > &space_dh, const Particles::ParticleHandler< dim, spacedim > &particle_handler, SparsityType &sparsity, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const ComponentMask &space_comps=ComponentMask())
Definition: utilities.cc:32
typename ActiveSelector::cell_iterator cell_iterator
Definition: dof_handler.h:466
particle_iterator end() const
static ::ExceptionBase & ExcInternalError()