using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
Number>::vectorized_value_type;
+ /**
+ * Collects the options which can be used to specify the behavior during
+ * reinitialization.
+ */
+ struct AdditionalData
+ {
+ /**
+ * Constructor which sets the default arguments.
+ */
+ AdditionalData(const bool use_global_weights = false)
+ : use_global_weights(use_global_weights)
+ {}
+
+ /**
+ * During initialization, assume that the Quadrature object contains
+ * global weights as, e.g., obtained by
+ * QSimplex::compute_affine_transformation().
+ */
+ bool use_global_weights;
+ };
+
/**
* Constructor.
*
* @param update_flags Specify the quantities to be computed by the mapping
* during the call of reinit(). These update flags are also handed to a
* FEEvaluation object if you construct it with this MappingInfo object.
+ *
+ * @param additional_data Additional data for the class to specify the
+ * behavior during reinitialization.
*/
- MappingInfo(const Mapping<dim> &mapping, const UpdateFlags update_flags);
+ MappingInfo(const Mapping<dim> & mapping,
+ const UpdateFlags update_flags,
+ const AdditionalData additional_data = AdditionalData());
/**
* Compute the mapping information for the incoming cell and unit
* Store the requested mapping data.
*/
void
- store_mapping_data(const unsigned int unit_points_index_offset,
- const unsigned int n_q_points,
- const unsigned int n_q_points_unvectorized,
- const MappingData &mapping_data);
+ store_mapping_data(const unsigned int unit_points_index_offset,
+ const unsigned int n_q_points,
+ const unsigned int n_q_points_unvectorized,
+ const MappingData & mapping_data,
+ const std::vector<double> &weights);
/**
* Compute the compressed cell index.
*/
UpdateFlags update_flags_mapping;
+ /**
+ * AdditionalData for this object.
+ */
+ const AdditionalData additional_data;
+
/**
* Stores the index offset into the arrays @p JxW_values, @p jacobians,
* @p inverse_jacobians and @p normal_vectors.
template <int dim, int spacedim, typename Number>
MappingInfo<dim, spacedim, Number>::MappingInfo(
- const Mapping<dim> &mapping,
- const UpdateFlags update_flags)
+ const Mapping<dim> & mapping,
+ const UpdateFlags update_flags,
+ const AdditionalData additional_data)
: mapping(&mapping)
, update_flags(update_flags)
, update_flags_mapping(update_default)
+ , additional_data(additional_data)
{
// translate update flags
if (update_flags & update_jacobians)
store_mapping_data(0,
n_q_points_data,
n_q_points_unvectorized[0],
- mapping_data);
+ mapping_data,
+ quadrature.get_weights());
state = State::single_cell;
is_reinitialized();
store_mapping_data(data_index_offsets[cell_index],
n_q_points_data,
n_q_points_unvectorized[cell_index],
- mapping_data);
+ mapping_data,
+ quadrature_vector[cell_index].get_weights());
if (do_cell_index_compression)
cell_index_to_compressed_cell_index[cell->active_cell_index()] =
const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
const unsigned int n_unfiltered_cells)
{
+ Assert(
+ additional_data.use_global_weights == false,
+ ExcMessage(
+ "There is no known use-case for AdditionalData::use_global_weights=true and reinit_surface()"));
+
do_cell_index_compression =
n_unfiltered_cells != numbers::invalid_unsigned_int;
store_mapping_data(data_index_offsets[cell_index],
n_q_points_data,
n_q_points_unvectorized[cell_index],
- mapping_data);
+ mapping_data,
+ quadrature_vector[cell_index].get_weights());
if (do_cell_index_compression)
cell_index_to_compressed_cell_index[cell->active_cell_index()] =
store_mapping_data(data_index_offsets[current_face_index],
n_q_points_data,
n_q_points_unvectorized[current_face_index],
- mapping_data);
+ mapping_data,
+ quadrature_on_face.get_weights());
}
if (do_cell_index_compression)
cell_index_to_compressed_cell_index[cell->active_cell_index()] =
const unsigned int unit_points_index_offset,
const unsigned int n_q_points,
const unsigned int n_q_points_unvectorized,
- const MappingInfo::MappingData &mapping_data)
+ const MappingInfo::MappingData &mapping_data,
+ const std::vector<double> & weights)
{
const unsigned int n_lanes =
dealii::internal::VectorizedArrayTrait<Number>::width();
inverse_jacobians[offset][s][d], v) =
mapping_data.inverse_jacobians[q * n_lanes + v][s][d];
if (update_flags_mapping & UpdateFlags::update_JxW_values)
- dealii::internal::VectorizedArrayTrait<Number>::get(
- JxW_values[offset], v) =
- mapping_data.JxW_values[q * n_lanes + v];
+ {
+ if (additional_data.use_global_weights)
+ {
+ dealii::internal::VectorizedArrayTrait<Number>::get(
+ JxW_values[offset], v) = weights[q * n_lanes + v];
+ }
+ else
+ {
+ dealii::internal::VectorizedArrayTrait<Number>::get(
+ JxW_values[offset], v) =
+ mapping_data.JxW_values[q * n_lanes + v];
+ }
+ }
if (update_flags_mapping & UpdateFlags::update_normal_vectors)
for (unsigned int s = 0; s < spacedim; ++s)
dealii::internal::VectorizedArrayTrait<Number>::get(
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+/*
+ * Test NonMatching::MappingInfo if reinit with quadratures which contain JxW as
+ * weights.
+ */
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_point_evaluation.h>
+
+#include <deal.II/non_matching/mapping_info.h>
+
+#include "../tests.h"
+
+void
+test()
+{
+ constexpr unsigned int dim = 2;
+ constexpr unsigned int spacedim = 2;
+ constexpr unsigned int degree = 1;
+ constexpr unsigned int n_components = 1;
+ using Number = double;
+
+ Triangulation<dim, spacedim> tria;
+ GridGenerator::subdivided_hyper_cube(tria, 4);
+
+ FE_Q<dim> fe(degree);
+ MappingQ<dim> mapping(degree);
+
+
+ deallog << "Check JxW faces..." << std::endl;
+ {
+ NonMatching::MappingInfo<dim>::AdditionalData additional_data;
+ additional_data.use_global_weights = true;
+ NonMatching::MappingInfo<dim> mapping_info(mapping,
+ update_JxW_values,
+ additional_data);
+
+ // 1) build vector of quadratures
+ std::vector<std::vector<Quadrature<dim - 1>>> quad_vec;
+ // prescribe JxW (there is no meaning in the actual values, they just have
+ // to stay the same when fetched with FEPointEvaluation)
+ double JxW = 1.0;
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ std::vector<Quadrature<dim - 1>> quad;
+ for (auto f : cell->face_indices())
+ {
+ dealii::QGauss<dim - 1> face_quadrature(degree + 1);
+ std::vector<double> weights(face_quadrature.get_weights().size());
+ for (auto &w : weights)
+ {
+ w = JxW;
+ JxW += 1.0;
+ }
+
+ quad.emplace_back(
+ Quadrature<dim - 1>(face_quadrature.get_points(), weights));
+ }
+ quad_vec.push_back(quad);
+ }
+
+ // 2) reinit mapping info
+ mapping_info.reinit_faces(tria.active_cell_iterators(), quad_vec);
+
+ FEPointEvaluation<n_components, dim, spacedim, Number> fe_point_eval(
+ mapping_info, fe);
+
+ // 3) print JxW
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ for (auto f : cell->face_indices())
+ {
+ fe_point_eval.reinit(cell->active_cell_index(), f);
+ for (unsigned int q : fe_point_eval.quadrature_point_indices())
+ deallog << fe_point_eval.JxW(q) << std::endl;
+ }
+ }
+ }
+
+ deallog << "\n\nCheck JxW cells..." << std::endl;
+ {
+ NonMatching::MappingInfo<dim>::AdditionalData additional_data;
+ additional_data.use_global_weights = true;
+ NonMatching::MappingInfo<dim> mapping_info(mapping,
+ update_JxW_values,
+ additional_data);
+
+
+ // 1) build vector of quadratures
+ std::vector<Quadrature<dim>> quad_vec;
+ // prescribe JxW (there is no meaning in the actual values, they just have
+ // to stay the same when fetched with FEPointEvaluation)
+ double JxW = 1.0;
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ dealii::QGauss<dim> cell_quadrature(degree + 1);
+ std::vector<double> weights(cell_quadrature.get_weights().size());
+ for (auto &w : weights)
+ {
+ w = JxW;
+ JxW += 1.0;
+ }
+
+ quad_vec.emplace_back(
+ Quadrature<dim>(cell_quadrature.get_points(), weights));
+ }
+
+ // 2) reinit mapping info
+ mapping_info.reinit_cells(tria.active_cell_iterators(), quad_vec);
+ FEPointEvaluation<n_components, dim, spacedim, Number> fe_point_eval(
+ mapping_info, fe);
+
+ // 3) print JxW
+ for (const auto &cell : tria.active_cell_iterators())
+ {
+ fe_point_eval.reinit(cell->active_cell_index());
+ for (unsigned int q : fe_point_eval.quadrature_point_indices())
+ deallog << fe_point_eval.JxW(q) << std::endl;
+ }
+ }
+
+
+ deallog << std::endl;
+}
+
+
+int
+main(int argc, char **argv)
+{
+ Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+ initlog();
+ test();
+}
projects / dealii.git / commitdiff