From: Johannes Heinz Date: Thu, 1 Jun 2023 15:44:45 +0000 (+0200) Subject: add NM::MappingInfo::AD to be able to reinit with global weights X-Git-Tag: v9.5.0-rc1~140^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5c32f7389a664d548a81b2dd7ba9ecc0c0ceb975;p=dealii.git add NM::MappingInfo::AD to be able to reinit with global weights --- diff --git a/doc/news/changes/minor/20230601Heinz b/doc/news/changes/minor/20230601Heinz new file mode 100644 index 0000000000..af0099c931 --- /dev/null +++ b/doc/news/changes/minor/20230601Heinz @@ -0,0 +1,3 @@ +Improved: Make it possible to fill JxW in NonMatching::MappingInfo from quadratures which contain JxW as weights. +
+(Johannes Heinz, 2023/06/01) diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h index 0f745240f3..85c7c512bd 100644 --- a/include/deal.II/non_matching/mapping_info.h +++ b/include/deal.II/non_matching/mapping_info.h @@ -207,6 +207,27 @@ namespace NonMatching 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. * @@ -215,8 +236,13 @@ namespace NonMatching * @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 &mapping, const UpdateFlags update_flags); + MappingInfo(const Mapping & mapping, + const UpdateFlags update_flags, + const AdditionalData additional_data = AdditionalData()); /** * Compute the mapping information for the incoming cell and unit @@ -448,10 +474,11 @@ namespace NonMatching * 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 &weights); /** * Compute the compressed cell index. @@ -523,6 +550,11 @@ namespace NonMatching */ 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. @@ -602,11 +634,13 @@ namespace NonMatching template MappingInfo::MappingInfo( - const Mapping &mapping, - const UpdateFlags update_flags) + const Mapping & 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) @@ -704,7 +738,8 @@ namespace NonMatching 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(); @@ -814,7 +849,8 @@ namespace NonMatching 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()] = @@ -837,6 +873,11 @@ namespace NonMatching const std::vector> &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; @@ -911,7 +952,8 @@ namespace NonMatching 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()] = @@ -1051,7 +1093,8 @@ namespace NonMatching 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()] = @@ -1219,7 +1262,8 @@ namespace NonMatching 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 & weights) { const unsigned int n_lanes = dealii::internal::VectorizedArrayTrait::width(); @@ -1244,9 +1288,19 @@ namespace NonMatching 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::get( - JxW_values[offset], v) = - mapping_data.JxW_values[q * n_lanes + v]; + { + if (additional_data.use_global_weights) + { + dealii::internal::VectorizedArrayTrait::get( + JxW_values[offset], v) = weights[q * n_lanes + v]; + } + else + { + dealii::internal::VectorizedArrayTrait::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::get( diff --git a/tests/non_matching/mapping_info_03.cc b/tests/non_matching/mapping_info_03.cc new file mode 100644 index 0000000000..cda12083aa --- /dev/null +++ b/tests/non_matching/mapping_info_03.cc @@ -0,0 +1,150 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include + +#include + +#include + +#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 tria; + GridGenerator::subdivided_hyper_cube(tria, 4); + + FE_Q fe(degree); + MappingQ mapping(degree); + + + deallog << "Check JxW faces..." << std::endl; + { + NonMatching::MappingInfo::AdditionalData additional_data; + additional_data.use_global_weights = true; + NonMatching::MappingInfo mapping_info(mapping, + update_JxW_values, + additional_data); + + // 1) build vector of quadratures + std::vector>> 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> quad; + for (auto f : cell->face_indices()) + { + dealii::QGauss face_quadrature(degree + 1); + std::vector weights(face_quadrature.get_weights().size()); + for (auto &w : weights) + { + w = JxW; + JxW += 1.0; + } + + quad.emplace_back( + Quadrature(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 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::AdditionalData additional_data; + additional_data.use_global_weights = true; + NonMatching::MappingInfo mapping_info(mapping, + update_JxW_values, + additional_data); + + + // 1) build vector of quadratures + std::vector> 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 cell_quadrature(degree + 1); + std::vector weights(cell_quadrature.get_weights().size()); + for (auto &w : weights) + { + w = JxW; + JxW += 1.0; + } + + quad_vec.emplace_back( + Quadrature(cell_quadrature.get_points(), weights)); + } + + // 2) reinit mapping info + mapping_info.reinit_cells(tria.active_cell_iterators(), quad_vec); + FEPointEvaluation 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(); +} diff --git a/tests/non_matching/mapping_info_03.output b/tests/non_matching/mapping_info_03.output new file mode 100644 index 0000000000..7c1ee20506 --- /dev/null +++ b/tests/non_matching/mapping_info_03.output @@ -0,0 +1,198 @@ + +DEAL::Check JxW faces... +DEAL::1.00000 +DEAL::2.00000 +DEAL::3.00000 +DEAL::4.00000 +DEAL::5.00000 +DEAL::6.00000 +DEAL::7.00000 +DEAL::8.00000 +DEAL::9.00000 +DEAL::10.0000 +DEAL::11.0000 +DEAL::12.0000 +DEAL::13.0000 +DEAL::14.0000 +DEAL::15.0000 +DEAL::16.0000 +DEAL::17.0000 +DEAL::18.0000 +DEAL::19.0000 +DEAL::20.0000 +DEAL::21.0000 +DEAL::22.0000 +DEAL::23.0000 +DEAL::24.0000 +DEAL::25.0000 +DEAL::26.0000 +DEAL::27.0000 +DEAL::28.0000 +DEAL::29.0000 +DEAL::30.0000 +DEAL::31.0000 +DEAL::32.0000 +DEAL::33.0000 +DEAL::34.0000 +DEAL::35.0000 +DEAL::36.0000 +DEAL::37.0000 +DEAL::38.0000 +DEAL::39.0000 +DEAL::40.0000 +DEAL::41.0000 +DEAL::42.0000 +DEAL::43.0000 +DEAL::44.0000 +DEAL::45.0000 +DEAL::46.0000 +DEAL::47.0000 +DEAL::48.0000 +DEAL::49.0000 +DEAL::50.0000 +DEAL::51.0000 +DEAL::52.0000 +DEAL::53.0000 +DEAL::54.0000 +DEAL::55.0000 +DEAL::56.0000 +DEAL::57.0000 +DEAL::58.0000 +DEAL::59.0000 +DEAL::60.0000 +DEAL::61.0000 +DEAL::62.0000 +DEAL::63.0000 +DEAL::64.0000 +DEAL::65.0000 +DEAL::66.0000 +DEAL::67.0000 +DEAL::68.0000 +DEAL::69.0000 +DEAL::70.0000 +DEAL::71.0000 +DEAL::72.0000 +DEAL::73.0000 +DEAL::74.0000 +DEAL::75.0000 +DEAL::76.0000 +DEAL::77.0000 +DEAL::78.0000 +DEAL::79.0000 +DEAL::80.0000 +DEAL::81.0000 +DEAL::82.0000 +DEAL::83.0000 +DEAL::84.0000 +DEAL::85.0000 +DEAL::86.0000 +DEAL::87.0000 +DEAL::88.0000 +DEAL::89.0000 +DEAL::90.0000 +DEAL::91.0000 +DEAL::92.0000 +DEAL::93.0000 +DEAL::94.0000 +DEAL::95.0000 +DEAL::96.0000 +DEAL::97.0000 +DEAL::98.0000 +DEAL::99.0000 +DEAL::100.000 +DEAL::101.000 +DEAL::102.000 +DEAL::103.000 +DEAL::104.000 +DEAL::105.000 +DEAL::106.000 +DEAL::107.000 +DEAL::108.000 +DEAL::109.000 +DEAL::110.000 +DEAL::111.000 +DEAL::112.000 +DEAL::113.000 +DEAL::114.000 +DEAL::115.000 +DEAL::116.000 +DEAL::117.000 +DEAL::118.000 +DEAL::119.000 +DEAL::120.000 +DEAL::121.000 +DEAL::122.000 +DEAL::123.000 +DEAL::124.000 +DEAL::125.000 +DEAL::126.000 +DEAL::127.000 +DEAL::128.000 +DEAL:: + +Check JxW cells... +DEAL::1.00000 +DEAL::2.00000 +DEAL::3.00000 +DEAL::4.00000 +DEAL::5.00000 +DEAL::6.00000 +DEAL::7.00000 +DEAL::8.00000 +DEAL::9.00000 +DEAL::10.0000 +DEAL::11.0000 +DEAL::12.0000 +DEAL::13.0000 +DEAL::14.0000 +DEAL::15.0000 +DEAL::16.0000 +DEAL::17.0000 +DEAL::18.0000 +DEAL::19.0000 +DEAL::20.0000 +DEAL::21.0000 +DEAL::22.0000 +DEAL::23.0000 +DEAL::24.0000 +DEAL::25.0000 +DEAL::26.0000 +DEAL::27.0000 +DEAL::28.0000 +DEAL::29.0000 +DEAL::30.0000 +DEAL::31.0000 +DEAL::32.0000 +DEAL::33.0000 +DEAL::34.0000 +DEAL::35.0000 +DEAL::36.0000 +DEAL::37.0000 +DEAL::38.0000 +DEAL::39.0000 +DEAL::40.0000 +DEAL::41.0000 +DEAL::42.0000 +DEAL::43.0000 +DEAL::44.0000 +DEAL::45.0000 +DEAL::46.0000 +DEAL::47.0000 +DEAL::48.0000 +DEAL::49.0000 +DEAL::50.0000 +DEAL::51.0000 +DEAL::52.0000 +DEAL::53.0000 +DEAL::54.0000 +DEAL::55.0000 +DEAL::56.0000 +DEAL::57.0000 +DEAL::58.0000 +DEAL::59.0000 +DEAL::60.0000 +DEAL::61.0000 +DEAL::62.0000 +DEAL::63.0000 +DEAL::64.0000 +DEAL::