From 13b0f1aab03c941913112a5f4d73f164227b3ed9 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Thu, 12 Sep 2024 14:29:52 +0200 Subject: [PATCH 1/4] Add new class MappingBox, implementing bbox mapping on general (hex or tet) cells --- examples/poisson.cc | 226 +++----- include/agglomeration_handler.h | 18 +- include/mapping_box.h | 385 +++++++++++++ source/CMakeLists.txt | 2 +- source/agglomeration_handler.cc | 221 +++---- source/mapping_box.cc | 988 ++++++++++++++++++++++++++++++++ 6 files changed, 1521 insertions(+), 319 deletions(-) create mode 100644 include/mapping_box.h create mode 100644 source/mapping_box.cc diff --git a/examples/poisson.cc b/examples/poisson.cc index 29f6e2e4..8266ec63 100644 --- a/examples/poisson.cc +++ b/examples/poisson.cc @@ -12,22 +12,19 @@ #include +#include #include #include #include #include -#include #include #include -#include #include #include #include -#include -#include #include #include @@ -35,7 +32,7 @@ #include #include - +#define HEX TRUE struct ConvergenceInfo { @@ -73,17 +70,6 @@ struct ConvergenceInfo -template -constexpr T -constexpr_pow(T num, unsigned int pow) -{ - return (pow >= sizeof(unsigned int) * 8) ? 0 : - pow == 0 ? 1 : - num * constexpr_pow(num, pow - 1); -} - - - enum class GridType { grid_generator, // hyper_cube or hyper_ball @@ -421,9 +407,14 @@ class Poisson output_results(); - Triangulation tria; - MappingQ mapping; - FE_DGQ dg_fe; + Triangulation tria; +#ifdef HEX + MappingQ1 mapping; + FE_DGQ dg_fe; +#else + MappingFE mapping; + FE_SimplexDGP dg_fe; +#endif std::unique_ptr> ah; AffineConstraints constraints; SparsityPattern sparsity; @@ -470,7 +461,12 @@ Poisson::Poisson(const GridType &grid_type, const unsigned int extraction_level, const unsigned int n_subdomains, const unsigned int fe_degree) - : mapping(1) + : +#ifdef HEX + mapping() +#else + mapping(FE_SimplexP{1}) +#endif , dg_fe(fe_degree) , grid_type(grid_type) , partitioner_type(partitioner_type) @@ -502,30 +498,42 @@ Poisson::make_grid() if constexpr (dim == 2) { grid_in.attach_triangulation(tria); +#ifdef HEX + std::ifstream gmsh_file( + "../../meshes/t3.msh"); // unstructured square made by triangles +#else std::ifstream gmsh_file( - "../../meshes/t3.msh"); // unstructured square [0,1]^2 + "../../meshes/square_simplex_coarser.msh"); // unstructured square + // made by triangles +#endif grid_in.read_msh(gmsh_file); - tria.refine_global(5); // 4 + tria.refine_global(2); // 4 } else if constexpr (dim == 3) { - // { - // grid_in.attach_triangulation(tria); - // std::ifstream gmsh_file("../../meshes/piston_3.inp"); // - // piston - // mesh grid_in.read_abaqus(gmsh_file); tria.refine_global(1); - // } grid_in.attach_triangulation(tria); - std::ifstream mesh_file( - "../../meshes/csf_brain_filled_centered_UCD.inp"); // piston mesh - grid_in.read_ucd(mesh_file); - tria.refine_global(1); +#ifdef HEX + std::ifstream filename("../../meshes/piston_3.inp"); // piston mesh + grid_in.read_abaqus(filename); +#else + std::ifstream filename( + "../../meshes/gray_level_image1.vtk"); // liver or brain domain + grid_in.read_vtk(filename); + +#endif } } else { +#ifdef HEX GridGenerator::hyper_cube(tria, 0., 1.); - tria.refine_global(8); + tria.refine_global(4); +#else + Triangulation tria_hex; + GridGenerator::hyper_cube(tria_hex, 0., 1.); + tria_hex.refine_global(4); + GridGenerator::convert_hypercube_to_simplex_mesh(tria_hex, tria); +#endif } std::cout << "Size of tria: " << tria.n_active_cells() << std::endl; cached_tria = std::make_unique>(tria, mapping); @@ -550,7 +558,7 @@ Poisson::make_grid() namespace bgi = boost::geometry::index; static constexpr unsigned int max_elem_per_node = - constexpr_pow(2, dim); // 2^dim + PolyUtils::constexpr_pow(2, dim); // 2^dim std::vector, typename Triangulation::active_cell_iterator>> boxes(tria.n_active_cells()); @@ -693,42 +701,6 @@ Poisson::setup_agglomeration() data_out.build_patches(mapping); data_out.write_vtu(output); } - - - - /* - #ifdef AGGLO_DEBUG - { - for (const auto &cell : ah->agglomeration_cell_iterators()) - { - std::cout << "Cell with idx: " << cell->active_cell_index() - << std::endl; - unsigned int n_agglomerated_faces_per_cell = ah->n_faces(cell); - std::cout << "Number of faces for the agglomeration: " - << n_agglomerated_faces_per_cell << std::endl; - for (unsigned int f = 0; f < n_agglomerated_faces_per_cell; ++f) - { - std::cout << "Agglomerated face with idx: " << f << std::endl; - const auto &[local_face_idx, neigh, local_face_idx_out, dummy] = - ah->get_agglomerated_connectivity()[{cell, f}]; - { - std::cout << "Face idx: " << local_face_idx << std::endl; - if (neigh.state() == IteratorState::valid) - { - std::cout << "Neighbor idx: " << - neigh->active_cell_index() - << std::endl; - } - - std::cout << "Face idx from outside: " << local_face_idx_out - << std::endl; - } - std::cout << std::endl; - } - } - } - #endif - */ } @@ -741,15 +713,24 @@ Poisson::assemble_system() solution.reinit(ah->n_dofs()); system_rhs.reinit(ah->n_dofs()); - const unsigned int quadrature_degree = 2 * dg_fe.get_degree() + 1; - const unsigned int face_quadrature_degree = 2 * dg_fe.get_degree() + 1; + const unsigned int quadrature_degree = dg_fe.get_degree() + 1; + const unsigned int face_quadrature_degree = dg_fe.get_degree() + 1; +#ifdef HEX ah->initialize_fe_values(QGauss(quadrature_degree), update_gradients | update_JxW_values | update_quadrature_points | update_JxW_values | update_values, QGauss(face_quadrature_degree)); +#else + ah->initialize_fe_values(QGaussSimplex(quadrature_degree), + update_gradients | update_JxW_values | + update_quadrature_points | update_JxW_values | + update_values, + QGaussSimplex(face_quadrature_degree)); +#endif const unsigned int dofs_per_cell = ah->n_dofs_per_cell(); + std::cout << "DoFs per cell: " << dofs_per_cell << std::endl; FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); Vector cell_rhs(dofs_per_cell); @@ -911,7 +892,7 @@ Poisson::assemble_system() ++i) { double d = (points0[i] - points1[i]).norm(); - Assert( + AssertThrow( d < 1e-15, ExcMessage( "Face qpoints at the interface do not match!")); @@ -930,11 +911,6 @@ Poisson::assemble_system() const auto &normals = fe_faces0.get_normal_vectors(); - // const double penalty = - // penalty_constant / - // PolyUtils::compute_h_orthogonal(f, - // polygon_boundary_vertices, - // normals[0]); const double penalty = penalty_constant / std::fabs(polytope->diameter()); @@ -1027,12 +1003,6 @@ Poisson::assemble_system() -void -output_double_number(double input, const std::string &text) -{ - std::cout << text << input << std::endl; -} - template void Poisson::solve() @@ -1048,27 +1018,6 @@ template void Poisson::output_results() { - { - const std::string filename = "agglomerated_Poisson.vtu"; - std::ofstream output(filename); - - DataOut data_out; - data_out.attach_dof_handler(ah->agglo_dh); - data_out.add_data_vector(solution, "u", DataOut::type_dof_data); - data_out.build_patches(/**(ah->euler_mapping)*/ mapping); - data_out.write_vtu(output); - } - { - const std::string filename = "agglomerated_Poisson_test.vtu"; - std::ofstream output(filename); - - DataOut data_out; - data_out.attach_dof_handler(ah->agglo_dh); - data_out.add_data_vector(solution, "u", DataOut::type_dof_data); - data_out.build_patches(*(ah->euler_mapping), 3); - data_out.write_vtu(output); - } - { std::string partitioner; if (partitioner_type == PartitionerType::metis) @@ -1078,13 +1027,16 @@ Poisson::output_results() else partitioner = "no_partitioning"; - const std::string filename = "interpolated_solution_" + partitioner + "_" + + const std::string filename = "interpolated_solution" + partitioner + "_" + std::to_string(n_subdomains) + ".vtu"; std::ofstream output(filename); DataOut data_out; Vector interpolated_solution; - PolyUtils::interpolate_to_fine_grid(*ah, interpolated_solution, solution); + PolyUtils::interpolate_to_fine_grid(*ah, + interpolated_solution, + solution, + true /*no_the_fly*/); data_out.attach_dof_handler(ah->output_dh); data_out.add_data_vector(interpolated_solution, "u", @@ -1120,53 +1072,6 @@ Poisson::output_results() semih1_err = errors[1]; std::cout << "Error (L2): " << l2_err << std::endl; std::cout << "Error (H1): " << semih1_err << std::endl; - -// Old version -#ifdef FALSE - { - // L2 error - Vector difference_per_cell(tria.n_active_cells()); - VectorTools::integrate_difference(mapping, - ah->output_dh, - interpolated_solution, - *analytical_solution, - difference_per_cell, - QGauss(dg_fe.degree), - VectorTools::L2_norm); - - // Plot the error per cell - - data_out.add_data_vector(difference_per_cell, - "error_per_cell", - DataOut::type_cell_data); - - const double L2_error = - VectorTools::compute_global_error(tria, - difference_per_cell, - VectorTools::L2_norm); - - std::cout << "L2 error:" << L2_error << std::endl; - } - - { - // H1 seminorm - Vector difference_per_cell(tria.n_active_cells()); - VectorTools::integrate_difference(mapping, - ah->output_dh, - interpolated_solution, - *analytical_solution, - difference_per_cell, - QGauss(dg_fe.degree + 1), - VectorTools::H1_seminorm); - - const double H1_seminorm = - VectorTools::compute_global_error(tria, - difference_per_cell, - VectorTools::H1_seminorm); - - std::cout << "H1 seminorm:" << H1_seminorm << std::endl; - } -#endif } } @@ -1196,7 +1101,14 @@ Poisson::run() { make_grid(); setup_agglomeration(); + auto start = std::chrono::high_resolution_clock::now(); assemble_system(); + auto stop = std::chrono::high_resolution_clock::now(); + auto duration = + std::chrono::duration_cast(stop - start); + + std::cout << "Time taken by assemble_system(): " << duration.count() / 1e6 + << " seconds" << std::endl; solve(); output_results(); } @@ -1210,7 +1122,11 @@ main() ConvergenceInfo convergence_info; std::cout << "Testing p-convergence" << std::endl; { - for (unsigned int fe_degree : {1, 2, 3, 4, 5, 6}) +#ifdef HEX + for (unsigned int fe_degree : {1, 2, 3, 4}) +#else + for (unsigned int fe_degree : {1, 2, 3}) +#endif { std::cout << "Fe degree: " << fe_degree << std::endl; Poisson<2> poisson_problem{GridType::unstructured, diff --git a/include/agglomeration_handler.h b/include/agglomeration_handler.h index c3c0748a..d4ab59f1 100644 --- a/include/agglomeration_handler.h +++ b/include/agglomeration_handler.h @@ -50,6 +50,7 @@ #include #include +#include #include #include @@ -467,11 +468,7 @@ class AgglomerationHandler : public Subscriptor */ DoFHandler output_dh; - - std::unique_ptr< - MappingFEField>> - euler_mapping; - + std::unique_ptr> box_mapping; /** * This function stores the information needed to identify which polytopes are @@ -535,6 +532,12 @@ class AgglomerationHandler : public Subscriptor initialize_agglomeration_data( const std::unique_ptr> &cache_tria); + void + update_agglomerate( + AgglomerationContainer &polytope, + const typename Triangulation::active_cell_iterator + &master_cell); + /** * Reinitialize the agglomeration data. */ @@ -562,8 +565,7 @@ class AgglomerationHandler : public Subscriptor * cells. This fills also the euler vector */ void - create_bounding_box(const AgglomerationContainer &polytope, - const types::global_cell_index master_idx); + create_bounding_box(const AgglomerationContainer &polytope); inline types::global_cell_index @@ -759,8 +761,6 @@ class AgglomerationHandler : public Subscriptor ObserverPointer> mapping; - std::unique_ptr> euler_fe; - std::unique_ptr> cached_tria; const MPI_Comm communicator; diff --git a/include/mapping_box.h b/include/mapping_box.h new file mode 100644 index 00000000..9dbfee85 --- /dev/null +++ b/include/mapping_box.h @@ -0,0 +1,385 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2001 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +#ifndef dealii_mapping_box_h +#define dealii_mapping_box_h + + +#include + +#include +#include + +#include + +#include + + +DEAL_II_NAMESPACE_OPEN + +/** + * @addtogroup mapping + * @{ + */ + +/** + * A class providing a mapping from the reference cell to cells that are + * axiparallel, i.e., that have the shape of rectangles (in 2d) or + * boxes (in 3d) with edges parallel to the coordinate directions. The + * class therefore provides functionality that is equivalent to what, + * for example, MappingQ would provide for such cells. However, knowledge + * of the shape of cells allows this class to be substantially more + * efficient. + * + * Specifically, the mapping is meant for cells for which the mapping from + * the reference to the real cell is a scaling along the coordinate + * directions: The transformation from reference coordinates $\hat {\mathbf + * x}$ to real coordinates $\mathbf x$ on each cell is of the form + * @f{align*}{ + * {\mathbf x}(\hat {\mathbf x}) + * = + * \begin{pmatrix} + * h_x & 0 \\ + * 0 & h_y + * \end{pmatrix} + * \hat{\mathbf x} + * + {\mathbf v}_0 + * @f} + * in 2d, and + * @f{align*}{ + * {\mathbf x}(\hat {\mathbf x}) + * = + * \begin{pmatrix} + * h_x & 0 & 0 \\ + * 0 & h_y & 0 \\ + * 0 & 0 & h_z + * \end{pmatrix} + * \hat{\mathbf x} + * + {\mathbf v}_0 + * @f} + * in 3d, where ${\mathbf v}_0$ is the bottom left vertex and $h_x,h_y,h_z$ + * are the extents of the cell along the axes. + * + * The class is intended for efficiency, and it does not do a whole lot of + * error checking. If you apply this mapping to a cell that does not conform + * to the requirements above, you will get strange results. + */ +template +class MappingBox : public Mapping +{ +public: + MappingBox(const std::vector> &local_boxes, + const std::map + &polytope_translator); + // for documentation, see the Mapping base class + virtual std::unique_ptr> + clone() const override; + + /** + * Return @p true because MappingBox preserves vertex + * locations. + */ + virtual bool + preserves_vertex_locations() const override; + + virtual bool + is_compatible_with(const ReferenceCell &reference_cell) const override; + + /** + * @name Mapping points between reference and real cells + * @{ + */ + + // for documentation, see the Mapping base class + virtual Point + transform_unit_to_real_cell( + const typename Triangulation::cell_iterator &cell, + const Point &p) const override; + + // for documentation, see the Mapping base class + virtual Point + transform_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const Point &p) const override; + + // for documentation, see the Mapping base class + virtual void + transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &real_points, + const ArrayView> &unit_points) const override; + + /** + * @} + */ + + /** + * @name Functions to transform tensors from reference to real coordinates + * @{ + */ + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + // for documentation, see the Mapping base class + virtual void + transform(const ArrayView> &input, + const MappingKind kind, + const typename Mapping::InternalDataBase &internal, + const ArrayView> &output) const override; + + /** + * @} + */ + + /** + * @name Interface with FEValues + * @{ + */ + + /** + * Storage for internal data of the mapping. See Mapping::InternalDataBase + * for an extensive description. + * + * This includes data that is computed once when the object is created (in + * get_data()) as well as data the class wants to store from between the + * call to fill_fe_values(), fill_fe_face_values(), or + * fill_fe_subface_values() until possible later calls from the finite + * element to functions such as transform(). The latter class of member + * variables are marked as 'mutable'. + */ + class InternalData : public Mapping::InternalDataBase + { + public: + /** + * Default constructor. + */ + InternalData() = default; + + /** + * Constructor that initializes the object with a quadrature. + */ + InternalData(const Quadrature &quadrature); + + // Documentation see Mapping::InternalDataBase. + virtual void + reinit(const UpdateFlags update_flags, + const Quadrature &quadrature) override; + + /** + * Return an estimate (in bytes) for the memory consumption of this object. + */ + virtual std::size_t + memory_consumption() const override; + + /** + * Extents of the last cell we have seen in the coordinate directions, + * i.e., hx, hy, hz. + */ + mutable Tensor<1, dim> cell_extents; + + /** + * Traslation term in F(\hat{x})=J\hat{x} + c. + */ + mutable Tensor<1, dim> traslation; + + /** + * Reciprocal of the extents of the last cell we have seen in the + * coordinate directions, i.e., hx, + * hy, hz. + */ + mutable Tensor<1, dim> inverse_cell_extents; + + /** + * The volume element + */ + mutable double volume_element; + + /** + * Location of quadrature points of faces or subfaces in 3d with all + * possible orientations. Can be accessed with the correct offset provided + * via QProjector::DataSetDescriptor. Not needed/used for cells. + */ + std::vector> quadrature_points; + }; + +private: + // documentation can be found in Mapping::requires_update_flags() + virtual UpdateFlags + requires_update_flags(const UpdateFlags update_flags) const override; + + // documentation can be found in Mapping::get_data() + virtual std::unique_ptr::InternalDataBase> + get_data(const UpdateFlags, const Quadrature &quadrature) const override; + + using Mapping::get_face_data; + + // documentation can be found in Mapping::get_subface_data() + virtual std::unique_ptr::InternalDataBase> + get_subface_data(const UpdateFlags flags, + const Quadrature &quadrature) const override; + + // documentation can be found in Mapping::fill_fe_values() + virtual CellSimilarity::Similarity + fill_fe_values( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const override; + + using Mapping::fill_fe_face_values; + + // documentation can be found in Mapping::fill_fe_subface_values() + virtual void + fill_fe_subface_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const override; + + // documentation can be found in Mapping::fill_fe_immersed_surface_values() + virtual void + fill_fe_immersed_surface_values( + const typename Triangulation::cell_iterator &cell, + const NonMatching::ImmersedSurfaceQuadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const override; + + /** + * @} + */ + + /** + * Update the cell_extents field of the incoming InternalData object with the + * size of the incoming cell. + */ + void + update_cell_extents( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const InternalData &data) const; + + /** + * Compute the quadrature points if the UpdateFlags of the incoming + * InternalData object say that they should be updated. + * + * Called from fill_fe_values. + */ + void + maybe_update_cell_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const InternalData &data, + const ArrayView> &unit_quadrature_points, + std::vector> &quadrature_points) const; + + /** + * Compute the normal vectors if the UpdateFlags of the incoming InternalData + * object say that they should be updated. + */ + void + maybe_update_normal_vectors( + const unsigned int face_no, + const InternalData &data, + std::vector> &normal_vectors) const; + + /** + * Since the Jacobian is constant for this mapping all derivatives of the + * Jacobian are identically zero. Fill these quantities with zeros if the + * corresponding update flags say that they should be updated. + */ + void + maybe_update_jacobian_derivatives( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const; + + + /** + * Compute the volume elements if the UpdateFlags of the incoming + * InternalData object say that they should be updated. + */ + void + maybe_update_volume_elements(const InternalData &data) const; + + /** + * Compute the Jacobians if the UpdateFlags of the incoming + * InternalData object say that they should be updated. + */ + void + maybe_update_jacobians( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const; + + /** + * Compute the inverse Jacobians if the UpdateFlags of the incoming + * InternalData object say that they should be updated. + */ + void + maybe_update_inverse_jacobians( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const; + + /** + * Vector of (local) bounding boxes + */ + std::vector> boxes; + + /** + * Map from global cell index to bounding box index + */ + std::map + polytope_translator; +}; + +/** @} */ + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index f45c6e7c..cac6a7e4 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -3,7 +3,7 @@ SET(_d2_build_types "Release;Debug") SET(Release_postfix "") SET(Debug_postfix ".g") -set(_files agglomeration_handler.cc multigrid_amg.cc) +set(_files agglomeration_handler.cc multigrid_amg.cc mapping_box.cc) FOREACH(_build_type ${_d2_build_types}) # Postfix to use everywhere diff --git a/source/agglomeration_handler.cc b/source/agglomeration_handler.cc index a412d103..f0221e3b 100644 --- a/source/agglomeration_handler.cc +++ b/source/agglomeration_handler.cc @@ -96,7 +96,7 @@ AgglomerationHandler::define_agglomerate( ++n_agglomerations; // an agglomeration has been performed, record it - create_bounding_box(cells, master_idx); // fill the vector of bboxes + create_bounding_box(cells); // fill the vector of bboxes // Finally, return a polygonal iterator to the polytope just constructed. return {cells[0], this}; @@ -202,10 +202,6 @@ AgglomerationHandler::initialize_agglomeration_data( mapping = &cache_tria->get_mapping(); agglo_dh.reinit(*tria); - FE_DGQ dummy_dg(1); - euler_fe = std::make_unique>(dummy_dg, spacedim); - euler_dh.reinit(*tria); - euler_dh.distribute_dofs(*euler_fe); if (const auto parallel_tria = dynamic_cast< const dealii::parallel::TriangulationBase *>(&*tria)) @@ -214,14 +210,10 @@ AgglomerationHandler::initialize_agglomeration_data( parallel_tria->global_active_cell_index_partitioner(); master_slave_relationships.reinit( cells_partitioner.lock()->locally_owned_range(), communicator); - - const IndexSet &local_eulerian_index_set = euler_dh.locally_owned_dofs(); - euler_vector.reinit(local_eulerian_index_set, communicator); } else { master_slave_relationships.reinit(tria->n_active_cells(), MPI_COMM_SELF); - euler_vector.reinit(euler_dh.n_dofs()); } polytope_cache.clear(); @@ -254,6 +246,9 @@ AgglomerationHandler::distribute_agglomerated_dofs( ExcNotImplemented( "Currently, this interface supports only DGQ and DGP bases.")); + box_mapping = std::make_unique>( + bboxes, + master2polygon); // construct bounding box mapping if (hybrid_mesh) { @@ -269,8 +264,9 @@ AgglomerationHandler::distribute_agglomerated_dofs( } - fe_collection.push_back(*fe); // master - fe_collection.push_back(FE_Nothing()); // slave + fe_collection.push_back(*fe); // master + fe_collection.push_back( + FE_Nothing(fe->reference_cell())); // slave initialize_hp_structure(); @@ -293,94 +289,18 @@ AgglomerationHandler::distribute_agglomerated_dofs( template void AgglomerationHandler::create_bounding_box( - const AgglomerationContainer &polytope, - const types::global_cell_index master_idx) + const AgglomerationContainer &polytope) { Assert(n_agglomerations > 0, ExcMessage("No agglomeration has been performed.")); Assert(dim > 1, ExcNotImplemented()); - std::vector dof_indices(euler_fe->dofs_per_cell); - std::vector> pts; // store all the vertices + std::vector> pts; // store all the vertices for (const auto &cell : polytope) for (const auto i : cell->vertex_indices()) pts.push_back(cell->vertex(i)); bboxes.emplace_back(pts); - - typename DoFHandler::cell_iterator polytope_dh(*polytope[0], - &euler_dh); - polytope_dh->get_dof_indices(dof_indices); - - - const auto &p0 = - bboxes[master2polygon.at(master_idx)].get_boundary_points().first; - const auto &p1 = - bboxes[master2polygon.at(master_idx)].get_boundary_points().second; - if constexpr (dim == 2) - { - euler_vector[dof_indices[0]] = p0[0]; - euler_vector[dof_indices[4]] = p0[1]; - // Lower right - euler_vector[dof_indices[1]] = p1[0]; - euler_vector[dof_indices[5]] = p0[1]; - // Upper left - euler_vector[dof_indices[2]] = p0[0]; - euler_vector[dof_indices[6]] = p1[1]; - // Upper right - euler_vector[dof_indices[3]] = p1[0]; - euler_vector[dof_indices[7]] = p1[1]; - } - else if constexpr (dim == 3) - { - // Lowers - - // left - euler_vector[dof_indices[0]] = p0[0]; - euler_vector[dof_indices[8]] = p0[1]; - euler_vector[dof_indices[16]] = p0[2]; - - // right - euler_vector[dof_indices[1]] = p1[0]; - euler_vector[dof_indices[9]] = p0[1]; - euler_vector[dof_indices[17]] = p0[2]; - - // left - euler_vector[dof_indices[2]] = p0[0]; - euler_vector[dof_indices[10]] = p1[1]; - euler_vector[dof_indices[18]] = p0[2]; - - // right - euler_vector[dof_indices[3]] = p1[0]; - euler_vector[dof_indices[11]] = p1[1]; - euler_vector[dof_indices[19]] = p0[2]; - - // Uppers - - // left - euler_vector[dof_indices[4]] = p0[0]; - euler_vector[dof_indices[12]] = p0[1]; - euler_vector[dof_indices[20]] = p1[2]; - - // right - euler_vector[dof_indices[5]] = p1[0]; - euler_vector[dof_indices[13]] = p0[1]; - euler_vector[dof_indices[21]] = p1[2]; - - // left - euler_vector[dof_indices[6]] = p0[0]; - euler_vector[dof_indices[14]] = p1[1]; - euler_vector[dof_indices[22]] = p1[2]; - - // right - euler_vector[dof_indices[7]] = p1[0]; - euler_vector[dof_indices[15]] = p1[1]; - euler_vector[dof_indices[23]] = p1[2]; - } - else - { - Assert(false, ExcNotImplemented()); - } } @@ -567,9 +487,9 @@ AgglomerationHandler::initialize_hp_structure() ExcMessage("No agglomeration has been performed.")); agglo_dh.distribute_dofs(fe_collection); - euler_mapping = std::make_unique< - MappingFEField>>( - euler_dh, euler_vector); + // euler_mapping = std::make_unique< + // MappingFEField>>( euler_dh, euler_vector); } @@ -579,9 +499,10 @@ const FEValues & AgglomerationHandler::reinit( const AgglomerationIterator &polytope) const { - Assert(euler_mapping, - ExcMessage("The mapping describing the physical element stemming from " - "agglomeration has not been set up.")); + // Assert(euler_mapping, + // ExcMessage("The mapping describing the physical element stemming + // from " + // "agglomeration has not been set up.")); const auto &deal_cell = polytope->as_dof_handler_iterator(agglo_dh); @@ -594,22 +515,9 @@ AgglomerationHandler::reinit( Quadrature agglo_quad = agglomerated_quadrature(agglo_cells, deal_cell); - const double bbox_measure = - bboxes[master2polygon.at(deal_cell->active_cell_index())].volume(); - - // Scale weights with the volume of the BBox. This way, the euler_mapping - // defining the BBOx doesn't alter them. - std::vector scaled_weights; - std::transform(agglo_quad.get_weights().begin(), - agglo_quad.get_weights().end(), - std::back_inserter(scaled_weights), - [&bbox_measure](const double w) { return w / bbox_measure; }); - - Quadrature scaled_quad(agglo_quad.get_points(), scaled_weights); - - agglomerated_scratch = std::make_unique(*euler_mapping, + agglomerated_scratch = std::make_unique(*box_mapping, fe_collection[0], - scaled_quad, + agglo_quad, agglomeration_flags); return agglomerated_scratch->reinit(deal_cell); } @@ -636,9 +544,10 @@ AgglomerationHandler::reinit( const AgglomerationIterator &polytope, const unsigned int face_index) const { - Assert(euler_mapping, - ExcMessage("The mapping describing the physical element stemming from " - "agglomeration has not been set up.")); + // Assert(euler_mapping, + // ExcMessage("The mapping describing the physical element stemming + // from " + // "agglomeration has not been set up.")); const auto &deal_cell = polytope->as_dof_handler_iterator(agglo_dh); Assert(is_master_cell(deal_cell), ExcMessage("This should be true.")); @@ -685,9 +594,9 @@ AgglomerationHandler::reinit_interface( !neigh_polytope->is_locally_owned()), ExcInternalError()); - const auto &cell = polytope_in->as_dof_handler_iterator(agglo_dh); - const auto &bbox = bboxes[master2polygon.at(cell->active_cell_index())]; - const double bbox_measure = bbox.volume(); + const auto &cell = polytope_in->as_dof_handler_iterator(agglo_dh); + const auto &bbox = bboxes[master2polygon.at(cell->active_cell_index())]; + // const double bbox_measure = bbox.volume(); const unsigned int neigh_rank = neigh_polytope->subdomain_id(); const CellId &neigh_id = neigh_polytope->id(); @@ -699,7 +608,7 @@ AgglomerationHandler::reinit_interface( const auto &JxWs = recv_jxws.at(neigh_rank).at({neigh_id, local_neigh}); - const std::vector> &normals = + std::vector> &normals = recv_normals.at(neigh_rank).at({neigh_id, local_neigh}); // Apply the necessary scalings due to the bbox. @@ -711,31 +620,35 @@ AgglomerationHandler::reinit_interface( return bbox.real_to_unit(p); }); - std::vector scale_factors(final_unit_q_points.size()); - std::vector scaled_weights(final_unit_q_points.size()); - std::vector> scaled_normals(final_unit_q_points.size()); - - // Since we received normal vectors from a neighbor, we have to swap the - // sign of the vector in order to have outward normals. + // std::vector scale_factors(final_unit_q_points.size()); + // std::vector scaled_weights(final_unit_q_points.size()); + // std::vector> scaled_normals(final_unit_q_points.size()); + + // Since we received normal vectors from a neighbor, we have to swap + // the + // // sign of the vector in order to have outward normals. + // for (unsigned int q = 0; q < final_unit_q_points.size(); ++q) + // { + // for (unsigned int direction = 0; direction < spacedim; ++direction) + // scaled_normals[q][direction] = + // normals[q][direction] * (bbox.side_length(direction)); + + // scaled_normals[q] *= -1; + + // scaled_weights[q] = + // (JxWs[q] * scaled_normals[q].norm()) / bbox_measure; + // scaled_normals[q] /= scaled_normals[q].norm(); + // } for (unsigned int q = 0; q < final_unit_q_points.size(); ++q) - { - for (unsigned int direction = 0; direction < spacedim; ++direction) - scaled_normals[q][direction] = - normals[q][direction] * (bbox.side_length(direction)); + normals[q] *= -1; - scaled_normals[q] *= -1; - - scaled_weights[q] = - (JxWs[q] * scaled_normals[q].norm()) / bbox_measure; - scaled_normals[q] /= scaled_normals[q].norm(); - } NonMatching::ImmersedSurfaceQuadrature surface_quad( - final_unit_q_points, scaled_weights, scaled_normals); + final_unit_q_points, JxWs, normals); agglomerated_isv = std::make_unique>( - *euler_mapping, *fe, surface_quad, agglomeration_face_flags); + *box_mapping, *fe, surface_quad, agglomeration_face_flags); agglomerated_isv->reinit(cell); @@ -919,7 +832,7 @@ namespace dealii const auto &bbox = handler.bboxes[handler.master2polygon.at(cell->active_cell_index())]; - const double bbox_measure = bbox.volume(); + // const double bbox_measure = bbox.volume(); CellId polytope_out_id; if (neigh_polytope.state() == IteratorState::valid) @@ -959,27 +872,27 @@ namespace dealii // change the normals, we multiply by the correct factors in // order to obtain the original normal after the call to // `reinit(cell)`. - std::vector scale_factors(q_points.size()); - std::vector scaled_weights(q_points.size()); - std::vector> scaled_normals(q_points.size()); - - for (unsigned int q = 0; q < q_points.size(); ++q) - { - for (unsigned int direction = 0; direction < spacedim; - ++direction) - scaled_normals[q][direction] = - normals[q][direction] * (bbox.side_length(direction)); - - scaled_weights[q] = - (JxWs[q] * scaled_normals[q].norm()) / bbox_measure; - scaled_normals[q] /= scaled_normals[q].norm(); - } + // std::vector scale_factors(q_points.size()); + // std::vector scaled_weights(q_points.size()); + // std::vector> scaled_normals(q_points.size()); + + // for (unsigned int q = 0; q < q_points.size(); ++q) + // { + // for (unsigned int direction = 0; direction < spacedim; + // ++direction) + // scaled_normals[q][direction] = + // normals[q][direction] * (bbox.side_length(direction)); + + // scaled_weights[q] = + // (JxWs[q] * scaled_normals[q].norm()) / bbox_measure; + // scaled_normals[q] /= scaled_normals[q].norm(); + // } for (const auto &point : unit_q_points) final_unit_q_points.push_back(point); - for (const auto &weight : scaled_weights) + for (const auto &weight : JxWs) final_weights.push_back(weight); - for (const auto &normal : scaled_normals) + for (const auto &normal : normals) final_normals.push_back(normal); } @@ -988,7 +901,7 @@ namespace dealii agglo_isv_ptr = std::make_unique>( - *(handler.euler_mapping), + *(handler.box_mapping), *(handler.fe), surface_quad, handler.agglomeration_face_flags); diff --git a/source/mapping_box.cc b/source/mapping_box.cc new file mode 100644 index 00000000..d0b06887 --- /dev/null +++ b/source/mapping_box.cc @@ -0,0 +1,988 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2001 - 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include +#include + +#include + +#include + +#include + + + +DEAL_II_NAMESPACE_OPEN + +DeclExceptionMsg( + ExcCellNotAssociatedWithBox, + "You are using MappingBox, but the incoming element is not associated with a" + "Bounding Box Cartesian."); + + + +/** + * Return whether the incoming element has a BoundingBox associated to it. + * Simplicial and quad-hex meshes are supported. + */ +template +bool +has_box(const CellType &cell, + const std::map + &translator) +{ + Assert((cell->reference_cell().is_hyper_cube() || + cell->reference_cell().is_simplex()), + ExcNotImplemented()); + Assert((translator.find(cell->active_cell_index()) != translator.cend()), + ExcCellNotAssociatedWithBox()); + + return true; +} + + + +template +MappingBox::MappingBox( + const std::vector> &input_boxes, + const std::map + &global_to_polytope) +{ + Assert(input_boxes.size() > 0, + ExcMessage("Invalid number of bounding boxes.")); + + // copy boxes and map + boxes.resize(input_boxes.size()); + for (unsigned int i = 0; i < input_boxes.size(); ++i) + boxes[i] = input_boxes[i]; + polytope_translator = global_to_polytope; +} + + + +template +MappingBox::InternalData::InternalData(const Quadrature &q) + : cell_extents(numbers::signaling_nan>()) + , traslation(numbers::signaling_nan>()) + , inverse_cell_extents(numbers::signaling_nan>()) + , volume_element(numbers::signaling_nan()) + , quadrature_points(q.get_points()) +{} + + + +template +void +MappingBox::InternalData::reinit(const UpdateFlags update_flags, + const Quadrature &) +{ + // store the flags in the internal data object so we can access them + // in fill_fe_*_values(). use the transitive hull of the required + // flags + this->update_each = update_flags; +} + + + +template +std::size_t +MappingBox::InternalData::memory_consumption() const +{ + return (Mapping::InternalDataBase::memory_consumption() + + MemoryConsumption::memory_consumption(cell_extents) + + MemoryConsumption::memory_consumption(traslation) + + MemoryConsumption::memory_consumption(inverse_cell_extents) + + MemoryConsumption::memory_consumption(volume_element)); +} + + + +template +bool +MappingBox::preserves_vertex_locations() const +{ + return true; +} + + + +template +bool +MappingBox::is_compatible_with( + const ReferenceCell &reference_cell) const +{ + Assert(dim == reference_cell.get_dimension(), + ExcMessage("The dimension of your mapping (" + + Utilities::to_string(dim) + + ") and the reference cell cell_type (" + + Utilities::to_string(reference_cell.get_dimension()) + + " ) do not agree.")); + + return reference_cell.is_hyper_cube() || reference_cell.is_simplex(); +} + + + +template +UpdateFlags +MappingBox::requires_update_flags(const UpdateFlags in) const +{ + // this mapping is pretty simple in that it can basically compute + // every piece of information wanted by FEValues without requiring + // computing any other quantities. boundary forms are one exception + // since they can be computed from the normal vectors without much + // further ado + UpdateFlags out = in; + if (out & update_boundary_forms) + out |= update_normal_vectors; + + return out; +} + + + +template +std::unique_ptr::InternalDataBase> +MappingBox::get_data(const UpdateFlags update_flags, + const Quadrature &q) const +{ + std::unique_ptr::InternalDataBase> data_ptr = + std::make_unique(); + data_ptr->reinit(requires_update_flags(update_flags), q); + + return data_ptr; +} + + + +template +std::unique_ptr::InternalDataBase> +MappingBox::get_subface_data( + const UpdateFlags update_flags, + const Quadrature &quadrature) const +{ + (void)update_flags; + (void)quadrature; + DEAL_II_NOT_IMPLEMENTED(); + return {}; +} + + + +template +void +MappingBox::update_cell_extents( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const InternalData &data) const +{ + // Compute start point and sizes along axes. The vertices to be looked at + // are 1, 2, 4 compared to the base vertex 0. + if (cell_similarity != CellSimilarity::translation) + { + const BoundingBox ¤t_box = + boxes[polytope_translator.at(cell->active_cell_index())]; + const std::pair, Point> &bdary_points = + current_box.get_boundary_points(); + + for (unsigned int d = 0; d < dim; ++d) + { + const double cell_extent_d = current_box.side_length(d); + data.cell_extents[d] = cell_extent_d; + + data.traslation[d] = + .5 * (bdary_points.first[d] + + bdary_points.second[d]); // midpoint of each interval + + Assert(cell_extent_d != 0., + ExcMessage("Cell does not appear to be Cartesian!")); + data.inverse_cell_extents[d] = 1. / cell_extent_d; + } + } +} + + + +namespace +{ + template + void + transform_quadrature_points( + const BoundingBox &box, + const ArrayView> &unit_quadrature_points, + std::vector> &quadrature_points) + { + for (unsigned int i = 0; i < quadrature_points.size(); ++i) + quadrature_points[i] = box.unit_to_real(unit_quadrature_points[i]); + } +} // namespace + + + +template +void +MappingBox::maybe_update_cell_quadrature_points( + const typename Triangulation::cell_iterator &cell, + const InternalData &data, + const ArrayView> &unit_quadrature_points, + std::vector> &quadrature_points) const +{ + if (data.update_each & update_quadrature_points) + transform_quadrature_points( + boxes[polytope_translator.at(cell->active_cell_index())], + unit_quadrature_points, + quadrature_points); +} + + + +template +void +MappingBox::maybe_update_normal_vectors( + const unsigned int face_no, + const InternalData &data, + std::vector> &normal_vectors) const +{ + // compute normal vectors. All normals on a face have the same value. + if (data.update_each & update_normal_vectors) + { + Assert(face_no < GeometryInfo::faces_per_cell, ExcInternalError()); + std::fill(normal_vectors.begin(), + normal_vectors.end(), + GeometryInfo::unit_normal_vector[face_no]); + } +} + + + +template +void +MappingBox::maybe_update_jacobian_derivatives( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + if (cell_similarity != CellSimilarity::translation) + { + if (data.update_each & update_jacobian_grads) + for (unsigned int i = 0; i < output_data.jacobian_grads.size(); ++i) + output_data.jacobian_grads[i] = DerivativeForm<2, dim, spacedim>(); + + if (data.update_each & update_jacobian_pushed_forward_grads) + for (unsigned int i = 0; + i < output_data.jacobian_pushed_forward_grads.size(); + ++i) + output_data.jacobian_pushed_forward_grads[i] = Tensor<3, spacedim>(); + + if (data.update_each & update_jacobian_2nd_derivatives) + for (unsigned int i = 0; + i < output_data.jacobian_2nd_derivatives.size(); + ++i) + output_data.jacobian_2nd_derivatives[i] = + DerivativeForm<3, dim, spacedim>(); + + if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives) + for (unsigned int i = 0; + i < output_data.jacobian_pushed_forward_2nd_derivatives.size(); + ++i) + output_data.jacobian_pushed_forward_2nd_derivatives[i] = + Tensor<4, spacedim>(); + + if (data.update_each & update_jacobian_3rd_derivatives) + for (unsigned int i = 0; + i < output_data.jacobian_3rd_derivatives.size(); + ++i) + output_data.jacobian_3rd_derivatives[i] = + DerivativeForm<4, dim, spacedim>(); + + if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives) + for (unsigned int i = 0; + i < output_data.jacobian_pushed_forward_3rd_derivatives.size(); + ++i) + output_data.jacobian_pushed_forward_3rd_derivatives[i] = + Tensor<5, spacedim>(); + } +} + + + +template +void +MappingBox::maybe_update_volume_elements( + const InternalData &data) const +{ + if (data.update_each & update_volume_elements) + { + double volume = data.cell_extents[0]; + for (unsigned int d = 1; d < dim; ++d) + volume *= data.cell_extents[d]; + data.volume_element = volume; + } +} + + + +template +void +MappingBox::maybe_update_jacobians( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + // "compute" Jacobian at the quadrature points, which are all the + // same + if (data.update_each & update_jacobians) + if (cell_similarity != CellSimilarity::translation) + for (unsigned int i = 0; i < output_data.jacobians.size(); ++i) + { + output_data.jacobians[i] = DerivativeForm<1, dim, spacedim>(); + for (unsigned int j = 0; j < dim; ++j) + output_data.jacobians[i][j][j] = data.cell_extents[j]; + } +} + + + +template +void +MappingBox::maybe_update_inverse_jacobians( + const InternalData &data, + const CellSimilarity::Similarity cell_similarity, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + // "compute" inverse Jacobian at the quadrature points, which are + // all the same + if (data.update_each & update_inverse_jacobians) + if (cell_similarity != CellSimilarity::translation) + for (unsigned int i = 0; i < output_data.inverse_jacobians.size(); ++i) + { + output_data.inverse_jacobians[i] = Tensor<2, dim>(); + for (unsigned int j = 0; j < dim; ++j) + output_data.inverse_jacobians[i][j][j] = + data.inverse_cell_extents[j]; + } +} + + + +template +CellSimilarity::Similarity +MappingBox::fill_fe_values( + const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); + + // convert data object to internal data for this class. fails with + // an exception if that is not possible + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + + update_cell_extents(cell, cell_similarity, data); + + maybe_update_cell_quadrature_points(cell, + data, + quadrature.get_points(), + output_data.quadrature_points); + + // compute Jacobian determinant. all values are equal and are the + // product of the local lengths in each coordinate direction + if (data.update_each & (update_JxW_values | update_volume_elements)) + if (cell_similarity != CellSimilarity::translation) + { + double J = data.cell_extents[0]; + for (unsigned int d = 1; d < dim; ++d) + J *= data.cell_extents[d]; + data.volume_element = J; + if (data.update_each & update_JxW_values) + for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i) + output_data.JxW_values[i] = quadrature.weight(i); + } + + + maybe_update_jacobians(data, cell_similarity, output_data); + maybe_update_jacobian_derivatives(data, cell_similarity, output_data); + maybe_update_inverse_jacobians(data, cell_similarity, output_data); + + return cell_similarity; +} + + + +template +void +MappingBox::fill_fe_subface_values( + const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const unsigned int subface_no, + const Quadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + (void)cell; + (void)face_no; + (void)subface_no; + (void)quadrature; + (void)internal_data; + (void)output_data; + DEAL_II_NOT_IMPLEMENTED(); +} + + + +template +void +MappingBox::fill_fe_immersed_surface_values( + const typename Triangulation::cell_iterator &cell, + const NonMatching::ImmersedSurfaceQuadrature &quadrature, + const typename Mapping::InternalDataBase &internal_data, + internal::FEValuesImplementation::MappingRelatedData + &output_data) const +{ + AssertDimension(dim, spacedim); + Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); + + // Convert data object to internal data for this class. Fails with an + // exception if that is not possible. + Assert(dynamic_cast(&internal_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(internal_data); + + + update_cell_extents(cell, CellSimilarity::none, data); + + maybe_update_cell_quadrature_points(cell, + data, + quadrature.get_points(), + output_data.quadrature_points); + + if (data.update_each & update_normal_vectors) + for (unsigned int i = 0; i < output_data.normal_vectors.size(); ++i) + output_data.normal_vectors[i] = quadrature.normal_vector(i); + + if (data.update_each & update_JxW_values) + for (unsigned int i = 0; i < output_data.JxW_values.size(); ++i) + output_data.JxW_values[i] = quadrature.weight(i); + + maybe_update_volume_elements(data); + maybe_update_jacobians(data, CellSimilarity::none, output_data); + maybe_update_jacobian_derivatives(data, CellSimilarity::none, output_data); + maybe_update_inverse_jacobians(data, CellSimilarity::none, output_data); +} + + + +template +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d = 0; d < dim; ++d) + output[i][d] = input[i][d] * data.inverse_cell_extents[d]; + return; + } + + case mapping_contravariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d = 0; d < dim; ++d) + output[i][d] = input[i][d] * data.cell_extents[d]; + return; + } + case mapping_piola: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d = 0; d < dim; ++d) + output[i][d] = + input[i][d] * data.cell_extents[d] / data.volume_element; + return; + } + default: + DEAL_II_NOT_IMPLEMENTED(); + } +} + + + +template +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = + input[i][d1][d2] * data.inverse_cell_extents[d2]; + return; + } + + case mapping_contravariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2]; + return; + } + + case mapping_covariant_gradient: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * + data.inverse_cell_extents[d2] * + data.inverse_cell_extents[d1]; + return; + } + + case mapping_contravariant_gradient: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1]; + return; + } + + case mapping_piola: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / + data.volume_element; + return; + } + + case mapping_piola_gradient: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1] / + data.volume_element; + return; + } + + default: + DEAL_II_NOT_IMPLEMENTED(); + } +} + + + +template +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = + input[i][d1][d2] * data.inverse_cell_extents[d2]; + return; + } + + case mapping_contravariant: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2]; + return; + } + + case mapping_covariant_gradient: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * + data.inverse_cell_extents[d2] * + data.inverse_cell_extents[d1]; + return; + } + + case mapping_contravariant_gradient: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1]; + return; + } + + case mapping_piola: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] / + data.volume_element; + return; + } + + case mapping_piola_gradient: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + + for (unsigned int i = 0; i < output.size(); ++i) + for (unsigned int d1 = 0; d1 < dim; ++d1) + for (unsigned int d2 = 0; d2 < dim; ++d2) + output[i][d1][d2] = input[i][d1][d2] * data.cell_extents[d2] * + data.inverse_cell_extents[d1] / + data.volume_element; + return; + } + + default: + DEAL_II_NOT_IMPLEMENTED(); + } +} + + + +template +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_covariant_gradient: + { + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int q = 0; q < output.size(); ++q) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = input[q][i][j][k] * + data.inverse_cell_extents[j] * + data.inverse_cell_extents[k]; + } + return; + } + default: + DEAL_II_NOT_IMPLEMENTED(); + } +} + + + +template +void +MappingBox::transform( + const ArrayView> &input, + const MappingKind mapping_kind, + const typename Mapping::InternalDataBase &mapping_data, + const ArrayView> &output) const +{ + AssertDimension(input.size(), output.size()); + Assert(dynamic_cast(&mapping_data) != nullptr, + ExcInternalError()); + const InternalData &data = static_cast(mapping_data); + + switch (mapping_kind) + { + case mapping_contravariant_hessian: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + + for (unsigned int q = 0; q < output.size(); ++q) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = input[q][i][j][k] * + data.cell_extents[i] * + data.inverse_cell_extents[j] * + data.inverse_cell_extents[k]; + } + return; + } + + case mapping_covariant_hessian: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + + for (unsigned int q = 0; q < output.size(); ++q) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = input[q][i][j][k] * + (data.inverse_cell_extents[i] * + data.inverse_cell_extents[j]) * + data.inverse_cell_extents[k]; + } + + return; + } + + case mapping_piola_hessian: + { + Assert(data.update_each & update_covariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_covariant_transformation")); + Assert(data.update_each & update_contravariant_transformation, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_contravariant_transformation")); + Assert(data.update_each & update_volume_elements, + typename FEValuesBase::ExcAccessToUninitializedField( + "update_volume_elements")); + + for (unsigned int q = 0; q < output.size(); ++q) + for (unsigned int i = 0; i < spacedim; ++i) + for (unsigned int j = 0; j < spacedim; ++j) + for (unsigned int k = 0; k < spacedim; ++k) + { + output[q][i][j][k] = + input[q][i][j][k] * + (data.cell_extents[i] / data.volume_element * + data.inverse_cell_extents[j]) * + data.inverse_cell_extents[k]; + } + + return; + } + + default: + DEAL_II_NOT_IMPLEMENTED(); + } +} + + + +template +Point +MappingBox::transform_unit_to_real_cell( + const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ + Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); + Assert(dim == spacedim, ExcNotImplemented()); + + return boxes[polytope_translator.at(cell->active_cell_index())].unit_to_real( + p); +} + + + +template +Point +MappingBox::transform_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const Point &p) const +{ + Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); + Assert(dim == spacedim, ExcNotImplemented()); + + return boxes[polytope_translator.at(cell->active_cell_index())].real_to_unit( + p); +} + + + +template +void +MappingBox::transform_points_real_to_unit_cell( + const typename Triangulation::cell_iterator &cell, + const ArrayView> &real_points, + const ArrayView> &unit_points) const +{ + Assert(has_box(cell, polytope_translator), ExcCellNotAssociatedWithBox()); + AssertDimension(real_points.size(), unit_points.size()); + + if (dim != spacedim) + DEAL_II_NOT_IMPLEMENTED(); + for (unsigned int i = 0; i < real_points.size(); ++i) + unit_points[i] = + boxes[polytope_translator.at(cell->active_cell_index())].real_to_unit( + real_points[i]); +} + + + +template +std::unique_ptr> +MappingBox::clone() const +{ + return std::make_unique>(*this); +} + + +//--------------------------------------------------------------------------- +// explicit instantiations +template class MappingBox<1>; +template class MappingBox<2>; +template class MappingBox<3>; + + +DEAL_II_NAMESPACE_CLOSE From c1e7de4592fdcfb26049dc5496ef6dc7d56fc32d Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Mon, 16 Sep 2024 19:36:16 +0200 Subject: [PATCH 2/4] Replace euler mapping with box mapping --- examples/minimal_SIP.cc | 682 ++++-------------- include/agglomeration_handler.h | 12 + .../coarse_operator_from_matrix_free.cc | 2 +- test/polydeal/distributed_injection_01.cc | 5 +- .../distributed_poisson_sanity_check_01.cc | 4 +- .../distributed_poisson_sanity_check_02.cc | 4 +- ...lly_distributed_poisson_sanity_check_01.cc | 4 +- ...lly_distributed_poisson_sanity_check_02.cc | 4 +- test/polydeal/minimal_SIP_Poisson.cc | 29 - test/polydeal/poisson_sanity_check_01.cc | 4 +- test/polydeal/poisson_sanity_check_02.cc | 2 +- test/polydeal/poisson_sanity_check_03.cc | 4 +- 12 files changed, 164 insertions(+), 592 deletions(-) diff --git a/examples/minimal_SIP.cc b/examples/minimal_SIP.cc index 1fc14f74..5f9d4827 100644 --- a/examples/minimal_SIP.cc +++ b/examples/minimal_SIP.cc @@ -1,80 +1,58 @@ -// ----------------------------------------------------------------------------- -// -// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later -// Copyright (C) XXXX - YYYY by the polyDEAL authors -// -// This file is part of the polyDEAL library. -// -// Detailed license information governing the source code -// can be found in LICENSE.md at the top level directory. -// -// ----------------------------------------------------------------------------- - - #include +#include #include +#include +#include +#include +#include +#include #include #include #include +#include #include #include #include #include +#include -constexpr double entry_tol = 1e-14; - -template -class RightHandSide : public Function -{ -public: - RightHandSide() - : Function() - {} - - virtual void - value_list(const std::vector> &points, - std::vector &values, - const unsigned int /*component*/) const override; -}; - - -template -void -RightHandSide::value_list(const std::vector> &points, - std::vector &values, - const unsigned int /*component*/) const -{ - for (unsigned int i = 0; i < values.size(); ++i) - values[i] = 8. * numbers::PI * numbers::PI * - std::sin(2. * numbers::PI * points[i][0]) * - std::sin(2. * numbers::PI * points[i][1]); -} +using namespace dealii; template -class MyFunction : public Function +class LinearFunction : public Function { public: + LinearFunction(const std::vector &coeffs) + { + Assert(coeffs.size() <= dim, ExcMessage("Wrong size!")); + coefficients.resize(coeffs.size()); + for (size_t i = 0; i < coeffs.size(); i++) + coefficients[i] = coeffs[i]; + } virtual double value(const Point &p, const unsigned int component = 0) const override; + std::vector coefficients; }; template double -MyFunction::value(const Point &p, const unsigned int) const +LinearFunction::value(const Point &p, const unsigned int) const { - return p[0]; + double value = 0.; + for (size_t i = 0; i < coefficients.size(); i++) + value += coefficients[i] * p[i]; + return value; } - template -class Poisson +class LaplaceOperator { private: void @@ -86,80 +64,56 @@ class Poisson void solve(); void - output_results(); + perform_sanity_check(); + - bool to_agglomerate; Triangulation tria; MappingQ mapping; FE_DGQ dg_fe; std::unique_ptr> ah; - // no hanging node in DG discretization, we define an AffineConstraints object - // so we can use the distribute_local_to_global() directly. - AffineConstraints constraints; - SparsityPattern sparsity; - DynamicSparsityPattern dsp; - SparseMatrix system_matrix; - Vector solution; - Vector system_rhs; - std::unique_ptr> cached_tria; - std::unique_ptr> rhs_function; + AffineConstraints constraints; + SparsityPattern sparsity; + DynamicSparsityPattern dsp; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; + std::unique_ptr> cached_tria; public: - constexpr Poisson(const bool = true); + LaplaceOperator(const unsigned int, const unsigned int fe_degree = 1); void run(); - inline const SparseMatrix & - get_matrix() - { - return system_matrix; - } - double penalty = 20.; + double penalty_constant = 10.; + unsigned int n_subdomains; }; template -constexpr Poisson::Poisson(const bool agglomerated_or_not) - : to_agglomerate(agglomerated_or_not) - , mapping(1) - , dg_fe(1) +LaplaceOperator::LaplaceOperator(const unsigned int n_subdomains, + const unsigned int fe_degree) + : mapping(1) + , dg_fe(fe_degree) + , n_subdomains(n_subdomains) {} template void -Poisson::make_grid() +LaplaceOperator::make_grid() { - GridGenerator::hyper_cube(tria, -1, 1); - if constexpr (dim == 2) - { - if (to_agglomerate) - { - tria.refine_global(3); - std::cout << "N elements " << tria.n_active_cells() << std::endl; - } - else - { - tria.refine_global(2); - std::cout << "N elements " << tria.n_active_cells() << std::endl; - } - } - else if constexpr (dim == 3) - { - if (to_agglomerate) - { - tria.refine_global(1); - } - } - std::cout << "Number of cells: " << tria.n_active_cells() << std::endl; - std::ofstream out("test_grid.vtk"); - GridOut grid_out; - grid_out.write_vtk(tria, out); - std::cout << "Grid written " << std::endl; - - cached_tria = std::make_unique>(tria, mapping); - rhs_function = std::make_unique>(); - + GridIn grid_in; + grid_in.attach_triangulation(tria); + std::ifstream gmsh_file("../../meshes/t3.msh"); // unstructured square [0,1]^2 + grid_in.read_msh(gmsh_file); + tria.refine_global(3); + cached_tria = std::make_unique>(tria, mapping); + // Partition the triangulation with graph partitioner. + + GridTools::partition_triangulation(n_subdomains, + tria, + SparsityTools::Partitioner::metis); + std::cout << "N subdomains: " << n_subdomains << std::endl; constraints.close(); } @@ -167,182 +121,18 @@ Poisson::make_grid() template void -Poisson::setup_agglomeration() - +LaplaceOperator::setup_agglomeration() { ah = std::make_unique>(*cached_tria); - if (to_agglomerate) - { - if constexpr (dim == 2) - { - std::vector idxs_to_be_agglomerated = {0, - 1, - 2, - 3}; - - std::vector::active_cell_iterator> - cells_to_be_agglomerated; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated, - cells_to_be_agglomerated); - - - std::vector idxs_to_be_agglomerated2 = {4, - 5, - 6, - 7}; - - std::vector::active_cell_iterator> - cells_to_be_agglomerated2; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated2, - cells_to_be_agglomerated2); - - - std::vector idxs_to_be_agglomerated3 = {8, - 9, - 10, - 11}; - std::vector::active_cell_iterator> - cells_to_be_agglomerated3; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated3, - cells_to_be_agglomerated3); - - - std::vector idxs_to_be_agglomerated4 = { - 12, 13, 14, 15}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated4; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated4, - cells_to_be_agglomerated4); - - std::vector idxs_to_be_agglomerated5 = { - 16, 17, 18, 19}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated5; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated5, - cells_to_be_agglomerated5); - - std::vector idxs_to_be_agglomerated6 = { - 20, 21, 22, 23}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated6; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated6, - cells_to_be_agglomerated6); - - std::vector idxs_to_be_agglomerated7 = { - 24, 25, 26, 27}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated7; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated7, - cells_to_be_agglomerated7); - - std::vector idxs_to_be_agglomerated8 = { - 28, 29, 30, 31}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated8; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated8, - cells_to_be_agglomerated8); - - std::vector idxs_to_be_agglomerated9 = { - 32, 33, 34, 35}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated9; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated9, - cells_to_be_agglomerated9); - - std::vector idxs_to_be_agglomerated10 = { - 36, 37, 38, 39}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated10; - PolyUtils::collect_cells_for_agglomeration( - tria, idxs_to_be_agglomerated10, cells_to_be_agglomerated10); - - std::vector idxs_to_be_agglomerated11 = { - 40, 41, 42, 43}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated11; - PolyUtils::collect_cells_for_agglomeration( - tria, idxs_to_be_agglomerated11, cells_to_be_agglomerated11); - - std::vector idxs_to_be_agglomerated12 = { - 44, 45, 46, 47}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated12; - PolyUtils::collect_cells_for_agglomeration( - tria, idxs_to_be_agglomerated12, cells_to_be_agglomerated12); - - std::vector idxs_to_be_agglomerated13 = { - 48, 49, 50, 51}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated13; - PolyUtils::collect_cells_for_agglomeration( - tria, idxs_to_be_agglomerated13, cells_to_be_agglomerated13); - - std::vector idxs_to_be_agglomerated14 = { - 52, 53, 54, 55}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated14; - PolyUtils::collect_cells_for_agglomeration( - tria, idxs_to_be_agglomerated14, cells_to_be_agglomerated14); - - std::vector idxs_to_be_agglomerated15 = { - 56, 57, 58, 59}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated15; - PolyUtils::collect_cells_for_agglomeration( - tria, idxs_to_be_agglomerated15, cells_to_be_agglomerated15); - - std::vector idxs_to_be_agglomerated16 = { - 60, 61, 62, 63}; //{36,37} - std::vector::active_cell_iterator> - cells_to_be_agglomerated16; - PolyUtils::collect_cells_for_agglomeration( - tria, idxs_to_be_agglomerated16, cells_to_be_agglomerated16); - - - - // Agglomerate the cells just stored - ah->define_agglomerate(cells_to_be_agglomerated); - ah->define_agglomerate(cells_to_be_agglomerated2); - ah->define_agglomerate(cells_to_be_agglomerated3); - ah->define_agglomerate(cells_to_be_agglomerated4); - ah->define_agglomerate(cells_to_be_agglomerated5); - ah->define_agglomerate(cells_to_be_agglomerated6); - ah->define_agglomerate(cells_to_be_agglomerated7); - ah->define_agglomerate(cells_to_be_agglomerated8); - ah->define_agglomerate(cells_to_be_agglomerated9); - ah->define_agglomerate(cells_to_be_agglomerated10); - ah->define_agglomerate(cells_to_be_agglomerated11); - ah->define_agglomerate(cells_to_be_agglomerated12); - ah->define_agglomerate(cells_to_be_agglomerated13); - ah->define_agglomerate(cells_to_be_agglomerated14); - ah->define_agglomerate(cells_to_be_agglomerated15); - ah->define_agglomerate(cells_to_be_agglomerated16); - } - else if constexpr (dim == 3) - { - std::vector idxs_to_be_agglomerated = { - 0, 1, 2, 3, 4, 5, 6, 7}; + std::vector::active_cell_iterator>> + cells_per_subdomain(n_subdomains); + for (const auto &cell : tria.active_cell_iterators()) + cells_per_subdomain[cell->subdomain_id()].push_back(cell); - std::vector::active_cell_iterator> - cells_to_be_agglomerated; - PolyUtils::collect_cells_for_agglomeration(tria, - idxs_to_be_agglomerated, - cells_to_be_agglomerated); + for (std::size_t i = 0; i < cells_per_subdomain.size(); ++i) + ah->define_agglomerate(cells_per_subdomain[i]); - // Agglomerate the cells just stored - ah->define_agglomerate(cells_to_be_agglomerated); - } - } ah->distribute_agglomerated_dofs(dg_fe); ah->create_agglomeration_sparsity_pattern(dsp); sparsity.copy_from(dsp); @@ -352,7 +142,7 @@ Poisson::setup_agglomeration() template void -Poisson::assemble_system() +LaplaceOperator::assemble_system() { system_matrix.reinit(sparsity); solution.reinit(ah->n_dofs()); @@ -379,14 +169,11 @@ Poisson::assemble_system() FullMatrix M22(dofs_per_cell, dofs_per_cell); std::vector local_dof_indices(dofs_per_cell); - std::vector local_dof_indices_bdary_cell( - dofs_per_cell); - std::vector local_dof_indices_neighbor( - dofs_per_cell); for (const auto &polytope : ah->polytope_iterators()) { - std::cout << "Polytope with idx: " << polytope->index() << std::endl; + // std::cout << "Cell with idx: " << cell->active_cell_index() << + // std::endl; cell_matrix = 0; cell_rhs = 0; const auto &agglo_values = ah->reinit(polytope); @@ -394,7 +181,6 @@ Poisson::assemble_system() const auto &q_points = agglo_values.get_quadrature_points(); const unsigned int n_qpoints = q_points.size(); std::vector rhs(n_qpoints); - rhs_function->value_list(q_points, rhs); for (unsigned int q_index : agglo_values.quadrature_point_indices()) { @@ -417,100 +203,75 @@ Poisson::assemble_system() // Face terms const unsigned int n_faces = polytope->n_faces(); - std::cout << "Number of (generalized) faces: " << n_faces << std::endl; + AssertThrow(n_faces > 0, + ExcMessage( + "Invalid element: at least 4 faces are required.")); + // const double agglo_measure = + // bboxes[cell->active_cell_index()].volume(); + // unsigned int n_jumps = 0; for (unsigned int f = 0; f < n_faces; ++f) { - const double hf = 1.; // cell->face(0)->measure(); -#ifdef AGGLO_DEBUG - std::cout << "Face measure: " << hf << std::endl; -#endif + const double current_element_diameter = + std::fabs(polytope->diameter()); + // const double penalty = + // penalty_constant * (1. / current_element_diameter); + if (polytope->at_boundary(f)) { - std::cout << "at boundary!" << std::endl; const auto &fe_face = ah->reinit(polytope, f); - { - std::cout << "I qpoints sul boundary sono: " << std::endl; - for (const auto &q : fe_face.get_quadrature_points()) - std::cout << q << std::endl; - } const unsigned int dofs_per_cell = fe_face.dofs_per_cell; - std::cout << "With dofs_per_cell =" << fe_face.dofs_per_cell - << std::endl; - + std::vector local_dof_indices_bdary_cell( + dofs_per_cell); // Get normal vectors seen from each agglomeration. - const auto &normals = fe_face.get_normal_vectors(); - cell_matrix = 0.; - for (unsigned int q_index : fe_face.quadrature_point_indices()) + // const auto &normals = fe_face.get_normal_vectors(); + cell_matrix = 0.; + for ([[maybe_unused]] unsigned int q_index : + fe_face.quadrature_point_indices()) { for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) { cell_matrix(i, j) += - 0. * - (-fe_face.shape_value(i, q_index) * - fe_face.shape_grad(j, q_index) * - normals[q_index] - - fe_face.shape_grad(i, q_index) * normals[q_index] * - fe_face.shape_value(j, q_index) + - (penalty / hf) * fe_face.shape_value(i, q_index) * - fe_face.shape_value(j, q_index)) * - fe_face.JxW(q_index); + 0.; // zero out, for this sanity check we need to + // neglect boundary contributions } } } // distribute DoFs polytope->get_dof_indices(local_dof_indices_bdary_cell); - system_matrix.add(local_dof_indices_bdary_cell, cell_matrix); + constraints.distribute_local_to_global( + cell_matrix, local_dof_indices_bdary_cell, system_matrix); } else { const auto &neigh_polytope = polytope->neighbor(f); -#ifdef AGGLO_DEBUG - std::cout << "Neighbor is " << neigh_polytope->index() - << std::endl; -#endif + + const double neigh_element_diameter = + std::fabs(neigh_polytope->diameter()); + const double penalty = + penalty_constant * + std::max(1. / current_element_diameter, + 1. / neigh_element_diameter); // Cinv still missing // This is necessary to loop over internal faces only once. if (polytope->index() < neigh_polytope->index()) { unsigned int nofn = polytope->neighbor_of_agglomerated_neighbor(f); -#ifdef AGGLO_DEBUG - std::cout << "Neighbor of neighbor is:" << nofn << std::endl; -#endif - - Assert(neigh_polytope->neighbor(nofn)->index() == - polytope->index(), - ExcMessage("Impossible.")); const auto &fe_faces = ah->reinit_interface(polytope, neigh_polytope, f, nofn); -#ifdef AGGLO_DEBUG - std::cout << "Reinited the interface:" << nofn << std::endl; -#endif const auto &fe_faces0 = fe_faces.first; const auto &fe_faces1 = fe_faces.second; -#ifdef AGGLO_DEBUG - std::cout << "Local from current: " << f << std::endl; - std::cout << "Local from neighbor: " << nofn << std::endl; - - std::cout << "Jump between " << polytope->index() << " and " - << neigh_polytope->index() << std::endl; - { - std::cout << "Qpoints from first polytope: " << std::endl; - for (const auto &q : fe_faces0.get_quadrature_points()) - std::cout << q << std::endl; - std::cout << "Qpoints from secondo polytope: " << std::endl; - for (const auto &q : fe_faces1.get_quadrature_points()) - std::cout << q << std::endl; - } -#endif + + std::vector + local_dof_indices_neighbor(dofs_per_cell); M11 = 0.; M12 = 0.; @@ -522,9 +283,6 @@ Poisson::assemble_system() for (unsigned int q_index : fe_faces0.quadrature_point_indices()) { -#ifdef AGGLO_DEBUG - std::cout << normals[q_index] << std::endl; -#endif for (unsigned int i = 0; i < dofs_per_cell; ++i) { for (unsigned int j = 0; j < dofs_per_cell; ++j) @@ -536,8 +294,7 @@ Poisson::assemble_system() 0.5 * fe_faces0.shape_grad(j, q_index) * normals[q_index] * fe_faces0.shape_value(i, q_index) + - (penalty / hf) * - fe_faces0.shape_value(i, q_index) * + (penalty)*fe_faces0.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces0.JxW(q_index); @@ -548,8 +305,7 @@ Poisson::assemble_system() 0.5 * fe_faces1.shape_grad(j, q_index) * normals[q_index] * fe_faces0.shape_value(i, q_index) - - (penalty / hf) * - fe_faces0.shape_value(i, q_index) * + (penalty)*fe_faces0.shape_value(i, q_index) * fe_faces1.shape_value(j, q_index)) * fe_faces1.JxW(q_index); @@ -561,8 +317,7 @@ Poisson::assemble_system() 0.5 * fe_faces0.shape_grad(j, q_index) * normals[q_index] * fe_faces1.shape_value(i, q_index) - - (penalty / hf) * - fe_faces1.shape_value(i, q_index) * + (penalty)*fe_faces1.shape_value(i, q_index) * fe_faces0.shape_value(j, q_index)) * fe_faces1.JxW(q_index); @@ -574,19 +329,18 @@ Poisson::assemble_system() 0.5 * fe_faces1.shape_grad(j, q_index) * normals[q_index] * fe_faces1.shape_value(i, q_index) + - (penalty / hf) * - fe_faces1.shape_value(i, q_index) * + (penalty)*fe_faces1.shape_value(i, q_index) * fe_faces1.shape_value(j, q_index)) * fe_faces1.JxW(q_index); } } } -#ifdef AGGLO_DEBUG - std::cout << "Neighbor is " << neigh_polytope->index() - << std::endl; -#endif // distribute DoFs accordingly + // std::cout << "Neighbor is " << + // neigh_cell->active_cell_index() + // << std::endl; + // Retrieve DoFs info from the cell iterator. neigh_polytope->get_dof_indices(local_dof_indices_neighbor); constraints.distribute_local_to_global(M11, @@ -611,223 +365,57 @@ Poisson::assemble_system() } -template -void -Poisson::solve() -{ - SparseDirectUMFPACK A_direct; - { - { - std::cout << "Sanity check:" << std::endl; - - // Sanity check: v(x,y)=x - Vector interpx(ah->get_dof_handler().n_dofs()); - MyFunction xfunction{}; - VectorTools::interpolate(*(ah->euler_mapping), - ah->get_dof_handler(), - xfunction, - interpx); - - double value = system_matrix.matrix_scalar_product(interpx, interpx); - std::cout << "Test = " << value << std::endl; - } - - if (to_agglomerate) - { - std::ofstream output_file("agglo_SIP_matrix.txt"); - system_matrix.print_formatted(output_file); - } - else - { - std::ofstream output_file("std_SIP_matrix.txt"); - system_matrix.print_formatted(output_file); - } - } - A_direct.initialize(system_matrix); - A_direct.vmult(solution, system_rhs); -} - - template void -Poisson::output_results() +LaplaceOperator::perform_sanity_check() { { - const std::string filename = "agglomerated_Poisson.vtu"; - std::ofstream output(filename); - - DataOut data_out; - data_out.attach_dof_handler(ah->get_dof_handler()); - data_out.add_data_vector(solution, "u", DataOut::type_dof_data); - data_out.build_patches(mapping); - data_out.write_vtu(output); - } - { - const std::string filename = "agglomerated_Poisson_test.vtu"; - std::ofstream output(filename); - - DataOut data_out; - data_out.attach_dof_handler(ah->get_dof_handler()); - data_out.add_data_vector(solution, "u", DataOut::type_dof_data); - data_out.build_patches(*(ah->euler_mapping)); - data_out.write_vtu(output); - } - { - const std::string filename = "interpolated_solution.vtu"; - std::ofstream output(filename); - - // Setup interpolation onto classical underlying DoFHandler - Vector interpolated_solution; - PolyUtils::interpolate_to_fine_grid(*ah, interpolated_solution, solution); - - DataOut data_out; - data_out.attach_dof_handler(ah->output_dh); - data_out.add_data_vector(interpolated_solution, - "u", - DataOut::type_dof_data); - data_out.build_patches(mapping); - data_out.write_vtu(output); - } - { - const std::string filename = "interpolated_solution.vtu"; - std::ofstream output(filename); - - DataOut data_out; - Vector interpolated_solution; - PolyUtils::interpolate_to_fine_grid(*ah, interpolated_solution, solution); - data_out.attach_dof_handler(ah->output_dh); - data_out.add_data_vector(interpolated_solution, - "u", - DataOut::type_dof_data); - data_out.build_patches(mapping); - data_out.write_vtu(output); + // Sanity check: v(x,y)=x + Vector interpx(ah->get_dof_handler().n_dofs()); + Vector interpxplusy(ah->get_dof_handler().n_dofs()); + Vector interpone(ah->get_dof_handler().n_dofs()); + LinearFunction xfunction{{1, 0}}; + LinearFunction xplusfunction{{1, 1}}; + VectorTools::interpolate(ah->get_agglomeration_mapping(), + ah->get_dof_handler(), + xfunction, + interpx); + VectorTools::interpolate(ah->get_agglomeration_mapping(), + ah->get_dof_handler(), + xplusfunction, + interpxplusy); + const double valuex = system_matrix.matrix_scalar_product(interpx, interpx); + const double valuexplusy = + system_matrix.matrix_scalar_product(interpxplusy, interpxplusy); + std::cout << "Test with f(x,y)=x:" << valuex << std::endl; + std::cout << "Test with f(x,y)=x+y:" << valuexplusy << std::endl; + + interpone = 1.; + double value_one = + system_matrix.matrix_scalar_product(interpone, interpone); + std::cout << "Test with 1: " << value_one << std::endl; } - - // { - // // Sanity check: v(x,y)=x - // Vector interpx(ah->agglo_dh.n_dofs()); - // MyFunction xfunction{}; - // VectorTools::interpolate(mapping, ah->agglo_dh, xfunction, interpx); - - // double value = system_matrix.matrix_scalar_product(interpx, interpx); - // std::cout << "Test = " << value << std::endl; - // } - - // Compute some quality metrics: - std::cout << "Computing quality metrics:" << std::endl; - auto metrics = PolyUtils::compute_quality_metrics(*ah); - - // uniformity factors - const auto &uf = std::get<0>(metrics); - std::cout << "Average Uniformity Factor: " - << std::accumulate(uf.begin(), uf.end(), 0.) / uf.size() - << std::endl; -#ifdef AGGLO_DEBUG - for (const auto &m : uf) - std::cout << m << std::endl; -#endif - - // circle ratios - const auto &cr = std::get<1>(metrics); - std::cout << "Average Circle Ratio: " - << std::accumulate(cr.begin(), cr.end(), 0.) / cr.size() - << std::endl; - - // box ratios - const auto &br = std::get<2>(metrics); - std::cout << "Average Box Ratio: " - << std::accumulate(br.begin(), br.end(), 0.) / br.size() - << std::endl; - - // overlap factor - const double of = std::get<3>(metrics); - std::cout << "Overlap factor: " << of << std::endl; - -#ifdef AGGLO_DEBUG - for (const auto &m : cr) - std::cout << m << std::endl; -#endif } template void -Poisson::run() +LaplaceOperator::run() { make_grid(); setup_agglomeration(); assemble_system(); - solve(); - output_results(); + perform_sanity_check(); } int -main(int argc, char *argv[]) +main() { - constexpr unsigned int dim = 2; - if (argc == 1) + for (const unsigned int n_agglomerates : {50, 100, 120, 300, 400, 800}) { - std::cout << "Running SIP with agglomerations." << std::endl; - Poisson poisson_problem; - poisson_problem.run(); - } - else if (argc == 2) - { - if (std::strcmp(argv[1], "agglo") == 0) - { - std::cout << "Running SIP with 4 agglomerated cells." << std::endl; - Poisson poisson_problem(true); // agglomerate - poisson_problem.run(); - } - else if (std::strcmp(argv[1], "standard") == 0) - { - std::cout << "Running SIP with standard grid." << std::endl; - Poisson poisson_problem(false); // standard SIPDG - poisson_problem.run(); - } - else if (std::strcmp(argv[1], "test_equality") == 0) - { - std::cout << "Running the two techniques and testing for equality" - << std::endl; - Poisson standard_problem(false); - standard_problem.run(); - const auto &standard_matrix = standard_problem.get_matrix(); - Poisson agglo_problem(true); - agglo_problem.run(); - const auto &agglo_matrix = agglo_problem.get_matrix(); - - std::cout - << "Comparing entries of the two matrices by subtracting them:" - << std::endl; - for (unsigned int i = 0; i < standard_matrix.m(); ++i) - for (unsigned int j = 0; j < standard_matrix.n(); ++j) - { - std::cout << std::fabs(standard_matrix.el(i, j) - - agglo_matrix.el(i, j)) - << "\t"; - Assert( - std::fabs(standard_matrix.el(i, j) - agglo_matrix.el(i, j)) < - entry_tol, - ExcMessage( - "Matrices are not equivalent up to machine precision. Code is buggy.")); - } - std::cout << std::endl; - } - else - { - Assert( - false, - ExcMessage( - "Please choose between `agglo` and `standard` way of assembling the matrix. If you want to test if the two approcahes are equal, use `test_equality` as command line argument.")); - } - } - else - { - Assert( - false, - ExcMessage( - "Invalid number of command line arguments, aborting the program.")); + LaplaceOperator<2> sanity_check{n_agglomerates}; + sanity_check.run(); } return 0; -} +} \ No newline at end of file diff --git a/include/agglomeration_handler.h b/include/agglomeration_handler.h index d4ab59f1..7532fcb9 100644 --- a/include/agglomeration_handler.h +++ b/include/agglomeration_handler.h @@ -300,6 +300,9 @@ class AgglomerationHandler : public Subscriptor inline const Mapping & get_mapping() const; + inline const MappingBox & + get_agglomeration_mapping() const; + inline const std::vector> & get_local_bboxes() const; @@ -897,6 +900,15 @@ AgglomerationHandler::get_mapping() const +template +inline const MappingBox & +AgglomerationHandler::get_agglomeration_mapping() const +{ + return *box_mapping; +} + + + template inline const Triangulation & AgglomerationHandler::get_triangulation() const diff --git a/test/polydeal/coarse_operator_from_matrix_free.cc b/test/polydeal/coarse_operator_from_matrix_free.cc index 0b00027c..d0651cd8 100644 --- a/test/polydeal/coarse_operator_from_matrix_free.cc +++ b/test/polydeal/coarse_operator_from_matrix_free.cc @@ -526,7 +526,7 @@ TestMGMatrix::agglomerate_and_compute_level_matrices() LinearAlgebra::distributed::Vector dst_coarse_op; dst_coarse_op.reinit(agglomeration_handler->agglo_dh.locally_owned_dofs(), comm); - VectorTools::interpolate(*(agglomeration_handler->euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), agglomeration_handler->agglo_dh, func, src_coarse); diff --git a/test/polydeal/distributed_injection_01.cc b/test/polydeal/distributed_injection_01.cc index 8a5022d3..5f1ffad0 100644 --- a/test/polydeal/distributed_injection_01.cc +++ b/test/polydeal/distributed_injection_01.cc @@ -354,7 +354,8 @@ DistributedHierarchyTest::setup_agglomerated_problem() ah_coarse->agglo_dh.locally_owned_dofs(); TrilinosWrappers::MPI::Vector interp_coarse(locally_owned_dofs_coarse, comm); - VectorTools::interpolate(*(ah_coarse->euler_mapping), + VectorTools::interpolate(ah_coarse->get_agglomeration_mapping(), + , ah_coarse->agglo_dh, func, interp_coarse); @@ -425,7 +426,7 @@ DistributedHierarchyTest::setup_agglomerated_problem() // Compute error: TrilinosWrappers::MPI::Vector interp_fine( dst.locally_owned_elements(), comm); // take parallel layout from dst - VectorTools::interpolate(*(ah_fine->euler_mapping), + VectorTools::interpolate(ah_fine->get_agglomeration_mapping(), ah_fine->agglo_dh, func, interp_fine); diff --git a/test/polydeal/distributed_poisson_sanity_check_01.cc b/test/polydeal/distributed_poisson_sanity_check_01.cc index 5de4fff7..f4353722 100644 --- a/test/polydeal/distributed_poisson_sanity_check_01.cc +++ b/test/polydeal/distributed_poisson_sanity_check_01.cc @@ -613,12 +613,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); diff --git a/test/polydeal/distributed_poisson_sanity_check_02.cc b/test/polydeal/distributed_poisson_sanity_check_02.cc index a8d874c9..17d87422 100644 --- a/test/polydeal/distributed_poisson_sanity_check_02.cc +++ b/test/polydeal/distributed_poisson_sanity_check_02.cc @@ -467,12 +467,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); diff --git a/test/polydeal/fully_distributed_poisson_sanity_check_01.cc b/test/polydeal/fully_distributed_poisson_sanity_check_01.cc index a8ed8b03..46ef8f3d 100644 --- a/test/polydeal/fully_distributed_poisson_sanity_check_01.cc +++ b/test/polydeal/fully_distributed_poisson_sanity_check_01.cc @@ -497,12 +497,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); diff --git a/test/polydeal/fully_distributed_poisson_sanity_check_02.cc b/test/polydeal/fully_distributed_poisson_sanity_check_02.cc index 24d530d7..d3ce5e20 100644 --- a/test/polydeal/fully_distributed_poisson_sanity_check_02.cc +++ b/test/polydeal/fully_distributed_poisson_sanity_check_02.cc @@ -483,12 +483,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(*(ah.euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); diff --git a/test/polydeal/minimal_SIP_Poisson.cc b/test/polydeal/minimal_SIP_Poisson.cc index 6d80e7c8..64bdc4af 100644 --- a/test/polydeal/minimal_SIP_Poisson.cc +++ b/test/polydeal/minimal_SIP_Poisson.cc @@ -71,8 +71,6 @@ class Poisson assemble_system(); void solve(); - void - output_results(); bool to_agglomerate; Triangulation tria; @@ -475,32 +473,6 @@ Poisson::solve() -template -void -Poisson::output_results() -{ - { - const std::string filename = "agglomerated_Poisson.vtu"; - std::ofstream output(filename); - - DataOut data_out; - data_out.attach_dof_handler(ah->get_dof_handler()); - data_out.add_data_vector(solution, "u", DataOut::type_dof_data); - data_out.build_patches(mapping); - data_out.write_vtu(output); - } - { - const std::string filename = "agglomerated_Poisson_test.vtu"; - std::ofstream output(filename); - - DataOut data_out; - data_out.attach_dof_handler(ah->get_dof_handler()); - data_out.add_data_vector(solution, "u", DataOut::type_dof_data); - data_out.build_patches(*(ah->euler_mapping)); - data_out.write_vtu(output); - } -} - template void Poisson::run() @@ -509,7 +481,6 @@ Poisson::run() setup_agglomeration(); assemble_system(); solve(); - output_results(); } template diff --git a/test/polydeal/poisson_sanity_check_01.cc b/test/polydeal/poisson_sanity_check_01.cc index 937f1855..860ae948 100644 --- a/test/polydeal/poisson_sanity_check_01.cc +++ b/test/polydeal/poisson_sanity_check_01.cc @@ -428,11 +428,11 @@ LaplaceOperator::perform_sanity_check() Vector interpone(ah->get_dof_handler().n_dofs()); LinearFunction xfunction{{1, 0}}; LinearFunction xplusfunction{{1, 1}}; - VectorTools::interpolate(*(ah->euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah->get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(*(ah->euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah->get_dof_handler(), xplusfunction, interpxplusy); diff --git a/test/polydeal/poisson_sanity_check_02.cc b/test/polydeal/poisson_sanity_check_02.cc index c73bd111..7291469f 100644 --- a/test/polydeal/poisson_sanity_check_02.cc +++ b/test/polydeal/poisson_sanity_check_02.cc @@ -389,7 +389,7 @@ LaplaceOperator::perform_sanity_check() for (const auto &func : functions) { Vector interp_vector(ah->get_dof_handler().n_dofs()); - VectorTools::interpolate(*(ah->euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah->get_dof_handler(), *func, interp_vector); diff --git a/test/polydeal/poisson_sanity_check_03.cc b/test/polydeal/poisson_sanity_check_03.cc index ac778d5e..669254d5 100644 --- a/test/polydeal/poisson_sanity_check_03.cc +++ b/test/polydeal/poisson_sanity_check_03.cc @@ -377,11 +377,11 @@ LaplaceOperator::perform_sanity_check() Vector interpone(ah->get_dof_handler().n_dofs()); LinearFunction xfunction{{1, 0}}; LinearFunction xplusfunction{{1, 1}}; - VectorTools::interpolate(*(ah->euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah->get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(*(ah->euler_mapping), + VectorTools::interpolate(ah->get_agglomeration_mapping(), ah->get_dof_handler(), xplusfunction, interpxplusy); From 36b3e38b127f6225d5b882a69d19936338d170fe Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Mon, 16 Sep 2024 20:08:14 +0200 Subject: [PATCH 3/4] Remove leftovers --- include/agglomeration_handler.h | 7 +------ test/polydeal/coarse_operator_from_matrix_free.cc | 2 +- test/polydeal/distributed_injection_01.cc | 1 - test/polydeal/distributed_poisson_sanity_check_01.cc | 4 ++-- test/polydeal/distributed_poisson_sanity_check_02.cc | 4 ++-- test/polydeal/fully_distributed_poisson_sanity_check_01.cc | 4 ++-- test/polydeal/fully_distributed_poisson_sanity_check_02.cc | 4 ++-- 7 files changed, 10 insertions(+), 16 deletions(-) diff --git a/include/agglomeration_handler.h b/include/agglomeration_handler.h index 7532fcb9..f57dd67c 100644 --- a/include/agglomeration_handler.h +++ b/include/agglomeration_handler.h @@ -361,7 +361,7 @@ class AgglomerationHandler : public Subscriptor void print_agglomeration(StreamType &out) { - for (const auto &cell : euler_dh.active_cell_iterators()) + for (const auto &cell : tria->active_cell_iterators()) out << "Cell with index: " << cell->active_cell_index() << " has associated value: " << master_slave_relationships[cell->global_active_cell_index()] @@ -774,11 +774,6 @@ class AgglomerationHandler : public Subscriptor hp::FECollection fe_collection; - /** - * DoFHandler for the physical space - */ - DoFHandler euler_dh; - /** * Eulerian vector describing the new cells obtained by the bounding boxes */ diff --git a/test/polydeal/coarse_operator_from_matrix_free.cc b/test/polydeal/coarse_operator_from_matrix_free.cc index d0651cd8..efc5f5fc 100644 --- a/test/polydeal/coarse_operator_from_matrix_free.cc +++ b/test/polydeal/coarse_operator_from_matrix_free.cc @@ -526,7 +526,7 @@ TestMGMatrix::agglomerate_and_compute_level_matrices() LinearAlgebra::distributed::Vector dst_coarse_op; dst_coarse_op.reinit(agglomeration_handler->agglo_dh.locally_owned_dofs(), comm); - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(agglomeration_handler->get_agglomeration_mapping(), agglomeration_handler->agglo_dh, func, src_coarse); diff --git a/test/polydeal/distributed_injection_01.cc b/test/polydeal/distributed_injection_01.cc index 5f1ffad0..e7bdeb71 100644 --- a/test/polydeal/distributed_injection_01.cc +++ b/test/polydeal/distributed_injection_01.cc @@ -355,7 +355,6 @@ DistributedHierarchyTest::setup_agglomerated_problem() TrilinosWrappers::MPI::Vector interp_coarse(locally_owned_dofs_coarse, comm); VectorTools::interpolate(ah_coarse->get_agglomeration_mapping(), - , ah_coarse->agglo_dh, func, interp_coarse); diff --git a/test/polydeal/distributed_poisson_sanity_check_01.cc b/test/polydeal/distributed_poisson_sanity_check_01.cc index f4353722..43e651bd 100644 --- a/test/polydeal/distributed_poisson_sanity_check_01.cc +++ b/test/polydeal/distributed_poisson_sanity_check_01.cc @@ -613,12 +613,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); diff --git a/test/polydeal/distributed_poisson_sanity_check_02.cc b/test/polydeal/distributed_poisson_sanity_check_02.cc index 17d87422..7765b7f1 100644 --- a/test/polydeal/distributed_poisson_sanity_check_02.cc +++ b/test/polydeal/distributed_poisson_sanity_check_02.cc @@ -467,12 +467,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); diff --git a/test/polydeal/fully_distributed_poisson_sanity_check_01.cc b/test/polydeal/fully_distributed_poisson_sanity_check_01.cc index 46ef8f3d..f1d38683 100644 --- a/test/polydeal/fully_distributed_poisson_sanity_check_01.cc +++ b/test/polydeal/fully_distributed_poisson_sanity_check_01.cc @@ -497,12 +497,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); diff --git a/test/polydeal/fully_distributed_poisson_sanity_check_02.cc b/test/polydeal/fully_distributed_poisson_sanity_check_02.cc index d3ce5e20..3561c789 100644 --- a/test/polydeal/fully_distributed_poisson_sanity_check_02.cc +++ b/test/polydeal/fully_distributed_poisson_sanity_check_02.cc @@ -483,12 +483,12 @@ main(int argc, char *argv[]) LinearFunction<2> xfunction{{1, 0}}; LinearFunction<2> xplusyfunction{{1, 1}}; - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xfunction, interpx); - VectorTools::interpolate(ah->get_agglomeration_mapping(), + VectorTools::interpolate(ah.get_agglomeration_mapping(), ah.get_dof_handler(), xplusyfunction, interpxplusy); From dbff07da829f6d68558a8ea940a5a4a0fab2a4a1 Mon Sep 17 00:00:00 2001 From: Marco Feder Date: Tue, 17 Sep 2024 11:26:21 +0200 Subject: [PATCH 4/4] Improve allocation of qpoints inside reinit() --- source/agglomeration_handler.cc | 59 ++++++++++----------------------- 1 file changed, 17 insertions(+), 42 deletions(-) diff --git a/source/agglomeration_handler.cc b/source/agglomeration_handler.cc index f0221e3b..cfa56b01 100644 --- a/source/agglomeration_handler.cc +++ b/source/agglomeration_handler.cc @@ -832,8 +832,6 @@ namespace dealii const auto &bbox = handler.bboxes[handler.master2polygon.at(cell->active_cell_index())]; - // const double bbox_measure = bbox.volume(); - CellId polytope_out_id; if (neigh_polytope.state() == IteratorState::valid) polytope_out_id = neigh_polytope->id(); @@ -847,53 +845,30 @@ namespace dealii std::vector final_weights; std::vector> final_normals; + const unsigned int expected_qpoints = + common_face.size() * handler.agglomeration_face_quad.size(); + final_unit_q_points.reserve(expected_qpoints); + final_weights.reserve(expected_qpoints); + final_normals.reserve(expected_qpoints); + for (const auto &[deal_cell, local_face_idx] : common_face) { handler.no_face_values->reinit(deal_cell, local_face_idx); - auto q_points = handler.no_face_values->get_quadrature_points(); - const auto &JxWs = handler.no_face_values->get_JxW_values(); + const auto &q_points = + handler.no_face_values->get_quadrature_points(); + const auto &JxWs = handler.no_face_values->get_JxW_values(); const auto &normals = handler.no_face_values->get_normal_vectors(); - std::vector> unit_q_points; - std::transform(q_points.begin(), - q_points.end(), - std::back_inserter(unit_q_points), - [&](const Point &p) { - return bbox.real_to_unit(p); - }); - - // Weights must be scaled with det(J)*|J^-t n| for each - // quadrature point. Use the fact that we are using a BBox, so - // the jacobi entries are the side_length in each direction. - // Unfortunately, normal vectors will be scaled internally by - // deal.II by using a covariant transformation. In order not to - // change the normals, we multiply by the correct factors in - // order to obtain the original normal after the call to - // `reinit(cell)`. - // std::vector scale_factors(q_points.size()); - // std::vector scaled_weights(q_points.size()); - // std::vector> scaled_normals(q_points.size()); - - // for (unsigned int q = 0; q < q_points.size(); ++q) - // { - // for (unsigned int direction = 0; direction < spacedim; - // ++direction) - // scaled_normals[q][direction] = - // normals[q][direction] * (bbox.side_length(direction)); - - // scaled_weights[q] = - // (JxWs[q] * scaled_normals[q].norm()) / bbox_measure; - // scaled_normals[q] /= scaled_normals[q].norm(); - // } - - for (const auto &point : unit_q_points) - final_unit_q_points.push_back(point); - for (const auto &weight : JxWs) - final_weights.push_back(weight); - for (const auto &normal : normals) - final_normals.push_back(normal); + const unsigned int n_qpoints_agglo = q_points.size(); + + for (unsigned int q = 0; q < n_qpoints_agglo; ++q) + { + final_unit_q_points.push_back(bbox.real_to_unit(q_points[q])); + final_weights.push_back(JxWs[q]); + final_normals.push_back(normals[q]); + } } NonMatching::ImmersedSurfaceQuadrature surface_quad(