10#include <deal.II/base/mpi.h>
11#include <deal.II/base/partitioner.h>
12#include <deal.II/base/vectorization.h>
13#include <deal.II/lac/la_parallel_vector.h>
38 std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
40 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
42 const unsigned int n_components);
55 template <
typename Number,
57 int simd_length = dealii::VectorizedArray<Number>::size()>
59 :
public dealii::LinearAlgebra::distributed::Vector<Number>
74 using ScalarVector = dealii::LinearAlgebra::distributed::Vector<Number>;
80 using ScalarVector::operator=;
89 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
107 unsigned int component)
const;
123 unsigned int component);
134 template <
typename Number2 = Number,
135 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
146 template <
typename Number2 = Number,
147 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
163 template <
typename Number2 = Number,
164 typename Tensor = dealii::Tensor<1, n_comp, Number2>>
172 template <
typename Number,
int n_comp,
int simd_length>
175 const std::shared_ptr<const dealii::Utilities::MPI::Partitioner>
182 auto vector_partitioner =
185 dealii::LinearAlgebra::distributed::Vector<Number>::reinit(
190 template <
typename Number,
int n_comp,
int simd_length>
192 ScalarVector &scalar_vector,
unsigned int component)
const
196 "Cannot extract from a vector with zero components."));
198 Assert(n_comp * scalar_vector.get_partitioner()->locally_owned_size() ==
199 this->get_partitioner()->locally_owned_size(),
200 dealii::ExcMessage(
"Called with a scalar_vector argument that has "
201 "incompatible local range."));
202 const auto local_size =
203 scalar_vector.get_partitioner()->locally_owned_size();
204 for (
unsigned int i = 0; i < local_size; ++i)
205 scalar_vector.local_element(i) =
206 this->local_element(i * n_comp + component);
207 scalar_vector.update_ghost_values();
211 template <
typename Number,
int n_comp,
int simd_length>
213 const ScalarVector &scalar_vector,
unsigned int component)
217 "Cannot insert into a vector with zero components."));
219 Assert(n_comp * scalar_vector.get_partitioner()->locally_owned_size() ==
220 this->get_partitioner()->locally_owned_size(),
221 dealii::ExcMessage(
"Called with a scalar_vector argument that has "
222 "incompatible local range."));
223 const auto local_size =
224 scalar_vector.get_partitioner()->locally_owned_size();
225 for (
unsigned int i = 0; i < local_size; ++i)
226 this->local_element(i * n_comp + component) =
227 scalar_vector.local_element(i);
232 template <
typename Number,
int n_comp,
int simd_length>
233 template <
typename Number2,
typename Tensor>
234 DEAL_II_ALWAYS_INLINE
inline Tensor
236 const unsigned int i)
const
238 static_assert(std::is_same<Number2, typename Tensor::value_type>::value,
239 "dummy type mismatch");
243 if constexpr (n_comp == 0)
246 if constexpr (std::is_same<Number, Number2>::value) {
249 for (
unsigned int d = 0; d < n_comp; ++d)
250 tensor[d] = this->local_element(i * n_comp + d);
252 }
else if constexpr (std::is_same<VectorizedArray, Number2>::value) {
255 std::array<
unsigned int, VectorizedArray::size()> indices;
256 for (
unsigned int k = 0; k < VectorizedArray::size(); ++k)
257 indices[k] = k * n_comp;
259 dealii::vectorized_load_and_transpose(
260 n_comp, this->begin() + i * n_comp, indices.data(), &tensor[0]);
271 template <
typename Number,
int n_comp,
int simd_length>
272 template <
typename Number2,
typename Tensor>
273 DEAL_II_ALWAYS_INLINE
inline Tensor
275 const unsigned int *js)
const
277 static_assert(std::is_same<Number2, typename Tensor::value_type>::value,
278 "dummy type mismatch");
282 if constexpr (n_comp == 0)
285 if constexpr (std::is_same<Number, Number2>::value) {
288 for (
unsigned int d = 0; d < n_comp; ++d)
289 tensor[d] = this->local_element(js[0] * n_comp + d);
291 }
else if constexpr (std::is_same<VectorizedArray, Number2>::value) {
294 std::array<
unsigned int, VectorizedArray::size()> indices;
295 for (
unsigned int k = 0; k < VectorizedArray::size(); ++k)
296 indices[k] = js[k] * n_comp;
298 dealii::vectorized_load_and_transpose(
299 n_comp, this->begin(), indices.data(), &tensor[0]);
310 template <
typename Number,
int n_comp,
int simd_length>
311 template <
typename Number2,
typename Tensor>
312 DEAL_II_ALWAYS_INLINE
inline void
314 const Tensor &tensor,
const unsigned int i)
316 static_assert(std::is_same<Number2, typename Tensor::value_type>::value,
317 "dummy type mismatch");
320 if constexpr (n_comp == 0)
323 if constexpr (std::is_same<Number, Number2>::value) {
326 for (
unsigned int d = 0; d < n_comp; ++d)
327 this->local_element(i * n_comp + d) = tensor[d];
329 }
else if constexpr (std::is_same<VectorizedArray, Number2>::value) {
332 std::array<
unsigned int, VectorizedArray::size()> indices;
333 for (
unsigned int k = 0; k < VectorizedArray::size(); ++k)
334 indices[k] = k * n_comp;
336 dealii::vectorized_transpose_and_store(
false,
340 this->begin() + i * n_comp);
void write_tensor(const Tensor &tensor, const unsigned int i)
void extract_component(ScalarVector &scalar_vector, unsigned int component) const
Tensor get_tensor(const unsigned int *js) const
void insert_component(const ScalarVector &scalar_vector, unsigned int component)
dealii::VectorizedArray< Number, simd_length > VectorizedArray
dealii::LinearAlgebra::distributed::Vector< Number > ScalarVector
void reinit_with_scalar_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &scalar_partitioner)
Tensor get_tensor(const unsigned int i) const
std::shared_ptr< const dealii::Utilities::MPI::Partitioner > create_vector_partitioner(const std::shared_ptr< const dealii::Utilities::MPI::Partitioner > &scalar_partitioner, const unsigned int n_components)
dealii::LinearAlgebra::distributed::Vector< Number > ScalarVector