From: Maximilian Bergbauer Date: Tue, 18 Jul 2023 08:04:54 +0000 (+0200) Subject: Implement paths for FE_Q_iso_Q1 in NonMatching X-Git-Tag: relicensing~596^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F15803%2Fhead;p=dealii.git Implement paths for FE_Q_iso_Q1 in NonMatching --- diff --git a/include/deal.II/non_matching/immersed_surface_quadrature.h b/include/deal.II/non_matching/immersed_surface_quadrature.h index 764a7ede6f..9c31556753 100644 --- a/include/deal.II/non_matching/immersed_surface_quadrature.h +++ b/include/deal.II/non_matching/immersed_surface_quadrature.h @@ -122,6 +122,12 @@ namespace NonMatching const std::vector & weights, const std::vector> &normals); + /** + * Clears weights, points and normals vectors. + */ + void + clear(); + /** * Extend the given formula by an additional quadrature point. * The point, weight and normal should be with respect to reference space, diff --git a/include/deal.II/non_matching/quadrature_generator.h b/include/deal.II/non_matching/quadrature_generator.h index e9f5e7b00d..4612d95f5f 100644 --- a/include/deal.II/non_matching/quadrature_generator.h +++ b/include/deal.II/non_matching/quadrature_generator.h @@ -201,6 +201,12 @@ namespace NonMatching const hp::QCollection<1> &quadratures1D, const AdditionalData & additional_data = AdditionalData()); + /** + * Clears the inside, outside and surface quadratures. + */ + void + clear_quadratures(); + /** * Construct immersed quadratures rules for the incoming level set * function over the BoundingBox. @@ -216,6 +222,14 @@ namespace NonMatching void generate(const Function &level_set, const BoundingBox &box); + /** + * Same as above but does not clear quadratures and appends it to the + * existing quadrature instead. + */ + void + generate_append(const Function & level_set, + const BoundingBox &box); + /** * Return the quadrature rule for the region * $\{x \in B : \psi(x) < 0 \}$ @@ -305,6 +319,12 @@ namespace NonMatching const hp::QCollection<1> &quadratures1D, const AdditionalData & additional_data = AdditionalData()); + /** + * Clears the inside, outside and surface quadratures. + */ + void + clear_quadratures(); + /** * Construct immersed quadratures rules for the incoming level set * function on a given face of the BoundingBox. @@ -322,6 +342,15 @@ namespace NonMatching const BoundingBox &box, const unsigned int face_index); + /** + * Same as above but does not clear quadratures and appends it to the + * existing quadrature instead. + */ + void + generate_append(const Function & level_set, + const BoundingBox &box, + const unsigned int face_index); + /** * Return the quadrature rule for the region * $\{x \in F : \psi(x) < 0 \}$ @@ -403,6 +432,12 @@ namespace NonMatching const hp::QCollection<1> &quadratures1D, const AdditionalData & additional_data = AdditionalData()); + /** + * Does nothing. Exists for compatibility reasons. + */ + void + clear_quadratures(); + /** * Construct immersed quadratures rules for the incoming level set * function on a given face of the BoundingBox. @@ -417,6 +452,15 @@ namespace NonMatching const BoundingBox<1> &box, const unsigned int face_index); + /** + * Same as above but does not clear quadratures and appends it to the + * existing quadrature instead. + */ + void + generate_append(const Function<1> & level_set, + const BoundingBox<1> &box, + const unsigned int face_index); + /** * @copydoc FaceQuadratureGenerator::get_inside_quadrature() */ @@ -518,6 +562,12 @@ namespace NonMatching generate(const typename Triangulation::active_cell_iterator &cell); private: + /** + * Construct immersed quadratures for FE_Q_iso_Q1. + */ + void + generate_fe_q_iso_q1(const BoundingBox &unit_box); + /** * Function that describes our level set function in reference space. */ @@ -575,6 +625,13 @@ namespace NonMatching const unsigned int face_index); private: + /** + * Construct immersed quadratures for FE_Q_iso_Q1. + */ + void + generate_fe_q_iso_q1(const BoundingBox &unit_box, + unsigned int face_index); + /** * Function that describes our level set function in reference space. */ @@ -701,6 +758,12 @@ namespace NonMatching */ ExtendableQuadrature(const Quadrature &quadrature); + /** + * Clears weights and points vectors. + */ + void + clear(); + /** * Add a point with an associated weight to the quadrature. */ @@ -760,6 +823,12 @@ namespace NonMatching ExtendableQuadrature & quadrature_by_definiteness(const Definiteness definiteness); + /** + * Clears all quadratures. + */ + void + clear(); + /** * Quadrature for the region $\{x \in B : \psi_i(x) < 0 \forall i \}$ of * the box, $B$. @@ -1396,6 +1465,26 @@ namespace NonMatching virtual void set_active_cell( const typename Triangulation::active_cell_iterator &cell) = 0; + + /** + * Set the dof values and the bounding box of the subcell the + * function should be evaluated on. Relevant for FE_Q_iso_Q1. + */ + virtual void + set_subcell(const std::vector &mask, + const BoundingBox & subcell_box) = 0; + + /** + * Returns flag indicating if the finite element is FE_Q_iso_Q1. + */ + virtual bool + is_fe_q_iso_q1() const = 0; + + /** + * Number of subdivisions of the FE_Q_iso_Q1 element. + */ + virtual unsigned int + n_subdivisions() const = 0; }; } // namespace DiscreteQuadratureGeneratorImplementation diff --git a/source/non_matching/immersed_surface_quadrature.cc b/source/non_matching/immersed_surface_quadrature.cc index 55b6deef15..7e71e2b834 100644 --- a/source/non_matching/immersed_surface_quadrature.cc +++ b/source/non_matching/immersed_surface_quadrature.cc @@ -38,6 +38,17 @@ namespace NonMatching + template + inline void + ImmersedSurfaceQuadrature::clear() + { + this->quadrature_points.clear(); + this->weights.clear(); + this->normals.clear(); + } + + + template void ImmersedSurfaceQuadrature::push_back( diff --git a/source/non_matching/mesh_classifier.cc b/source/non_matching/mesh_classifier.cc index 743342bad1..6134b58b9c 100644 --- a/source/non_matching/mesh_classifier.cc +++ b/source/non_matching/mesh_classifier.cc @@ -17,6 +17,7 @@ #include +#include "deal.II/fe/fe_q_iso_q1.h" #include #include #include @@ -410,6 +411,24 @@ namespace NonMatching face_index, local_levelset_values); + const FiniteElement &fe = + level_set_description->get_fe_collection()[fe_index]; + + const FE_Q_iso_Q1 *fe_q_iso_q1 = + dynamic_cast *>(&fe); + + const FE_Poly *fe_poly = dynamic_cast *>(&fe); + + const bool is_linear = fe_q_iso_q1 != nullptr || + (fe_poly != nullptr && fe_poly->get_degree() == 1); + + // shortcut for linear elements + if (is_linear) + { + return internal::MeshClassifierImplementation::location_from_dof_signs( + local_levelset_values); + } + lagrange_to_bernstein_face[fe_index][face_index].solve( local_levelset_values); @@ -466,7 +485,8 @@ namespace NonMatching for (unsigned int i = 0; i < fe_collection.size(); i++) { const FiniteElement &element = fe_collection[i]; - const FE_Q *fe_q = dynamic_cast *>(&element); + const FE_Q_Base * fe_q = + dynamic_cast *>(&element); Assert(fe_q != nullptr, ExcNotImplemented()); const FE_Bernstein fe_bernstein(fe_q->get_degree()); diff --git a/source/non_matching/quadrature_generator.cc b/source/non_matching/quadrature_generator.cc index df5e9bb56c..6e908042fe 100644 --- a/source/non_matching/quadrature_generator.cc +++ b/source/non_matching/quadrature_generator.cc @@ -15,6 +15,9 @@ #include +#include "deal.II/fe/fe_q.h" +#include "deal.II/fe/fe_q_iso_q1.h" + #include #include @@ -641,6 +644,16 @@ namespace NonMatching + template + inline void + ExtendableQuadrature::clear() + { + this->quadrature_points.clear(); + this->weights.clear(); + } + + + template void ExtendableQuadrature::push_back(const Point &point, @@ -670,6 +683,18 @@ namespace NonMatching + template + void + QPartitioning::clear() + { + positive.clear(); + negative.clear(); + indefinite.clear(); + surface.clear(); + } + + + /** * Takes a (dim-1)-dimensional point from the cross-section (orthogonal * to direction) of the box. Creates the two dim-dimensional points, which @@ -870,7 +895,7 @@ namespace NonMatching void QGeneratorBase::clear_quadratures() { - q_partitioning = QPartitioning(); + q_partitioning.clear(); } @@ -1328,6 +1353,25 @@ namespace NonMatching set_active_cell(const typename Triangulation::active_cell_iterator &cell) override; + /** + * @copydoc CellWiseFunction::set_subcell() + */ + void + set_subcell(const std::vector &mask, + const BoundingBox & subcell_box) override; + + /** + * @copydoc CellWiseFunction::is_fe_q_iso_q1() + */ + bool + is_fe_q_iso_q1() const override; + + /** + * @copydoc CellWiseFunction::n_subdivisions() + */ + unsigned int + n_subdivisions() const override; + /** * @copydoc Function::value() * @@ -1386,16 +1430,33 @@ namespace NonMatching SmartPointer> element; /** - * DOF-indices of the cell in the last call to set_active_cell() + * DOF-indices of the cell in the last call to set_active_cell(). */ std::vector local_dof_indices; /** * Local solution values of the cell in the last call to - * set_active_cell() + * set_active_cell(). */ std::vector local_dof_values; + /** + * Local solution values of the subcell after the last call to + * set_subcell(). + */ + std::vector local_dof_values_subcell; + + /** + * Bounding box of the subcell after the last call to set_subcell(). + */ + BoundingBox subcell_box; + + /** + * Number of subdivisions per line of the FE_Q_iso_Q1 element. Set to + * numbers::invalid_unsigned_int for other elements. + */ + unsigned int n_subdivisions_per_line; + /** * Description of the 1d polynomial basis for tensor product elements * used for the fast path of this class using tensor product @@ -1422,6 +1483,7 @@ namespace NonMatching const VectorType & dof_values) : dof_handler(&dof_handler) , global_dof_values(&dof_values) + , n_subdivisions_per_line(numbers::invalid_unsigned_int) { Assert(dof_handler.n_dofs() == dof_values.size(), ExcDimensionMismatch(dof_handler.n_dofs(), dof_values.size())); @@ -1443,6 +1505,8 @@ namespace NonMatching const auto dof_handler_cell = cell->as_dof_handler_iterator(*dof_handler); + const FE_Q fe_q1(1); + // Save the element and the local dof values, since this is what we need // to evaluate the function. @@ -1454,9 +1518,22 @@ namespace NonMatching poly.clear(); element = &dof_handler_cell->get_fe(); - if (element->n_base_elements() == 1 && - dealii::internal::FEPointEvaluation::is_fast_path_supported( - *element, 0)) + const FiniteElement *fe = &dof_handler_cell->get_fe(); + + if (const FE_Q_iso_Q1 *fe_q_iso_q1 = + dynamic_cast *>( + &dof_handler_cell->get_fe())) + { + this->n_subdivisions_per_line = fe_q_iso_q1->get_degree(); + fe = &fe_q1; + local_dof_values_subcell.resize(fe_q1.n_dofs_per_cell()); + } + else + this->n_subdivisions_per_line = numbers::invalid_unsigned_int; + + if (fe->n_base_elements() == 1 && + dealii::internal::FEPointEvaluation::is_fast_path_supported(*fe, + 0)) { dealii::internal::MatrixFreeFunctions::ShapeInfo shape_info; @@ -1465,7 +1542,7 @@ namespace NonMatching renumber = shape_info.lexicographic_numbering; poly = dealii::internal::FEPointEvaluation::get_polynomial_space( - element->base_element(0)); + fe->base_element(0)); polynomials_are_hat_functions = (poly.size() == 2 && poly[0].value(0.) == 1. && @@ -1489,6 +1566,38 @@ namespace NonMatching + template + void + RefSpaceFEFieldFunction::set_subcell( + const std::vector &mask, + const BoundingBox & subcell_box) + { + for (unsigned int i = 0; i < local_dof_values_subcell.size(); ++i) + local_dof_values_subcell[i] = local_dof_values[renumber[mask[i]]]; + + this->subcell_box = subcell_box; + } + + + + template + bool + RefSpaceFEFieldFunction::is_fe_q_iso_q1() const + { + return n_subdivisions_per_line != numbers::invalid_unsigned_int; + } + + + + template + unsigned int + RefSpaceFEFieldFunction::n_subdivisions() const + { + return n_subdivisions_per_line; + } + + + template bool RefSpaceFEFieldFunction::cell_is_set() const @@ -1514,10 +1623,11 @@ namespace NonMatching // TODO: this could be extended to a component that is not zero return dealii::internal::evaluate_tensor_product_value( poly, - local_dof_values, - point, + this->is_fe_q_iso_q1() ? local_dof_values_subcell : + local_dof_values, + this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) : point, polynomials_are_hat_functions, - renumber); + this->is_fe_q_iso_q1() ? std::vector() : renumber); } else { @@ -1546,10 +1656,13 @@ namespace NonMatching // TODO: this could be extended to a component that is not zero return dealii::internal::evaluate_tensor_product_value_and_gradient( poly, - local_dof_values, - point, + this->is_fe_q_iso_q1() ? local_dof_values_subcell : + local_dof_values, + this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) : + point, polynomials_are_hat_functions, - renumber) + this->is_fe_q_iso_q1() ? std::vector() : + renumber) .second; } else @@ -1578,7 +1691,11 @@ namespace NonMatching { // TODO: this could be extended to a component that is not zero return dealii::internal::evaluate_tensor_product_hessian( - poly, local_dof_values, point, renumber); + poly, + this->is_fe_q_iso_q1() ? local_dof_values_subcell : + local_dof_values, + this->is_fe_q_iso_q1() ? subcell_box.real_to_unit(point) : point, + this->is_fe_q_iso_q1() ? std::vector() : renumber); } else { @@ -1591,6 +1708,57 @@ namespace NonMatching return symmetrize(hessian); } } + + + + template + BoundingBox + create_subcell_box(const BoundingBox & unit_box, + const std::array &subcell_indices, + const unsigned int n_subdivisions) + { + std::pair, Point> corners; + for (unsigned int d = 0; d < dim; ++d) + { + AssertIndexRange(subcell_indices[d], n_subdivisions + 1); + + const double split_box_side_length = + (1. / n_subdivisions) * unit_box.side_length(d); + corners.first[d] = subcell_indices[d] * split_box_side_length; + corners.second[d] = corners.first[d] + split_box_side_length; + } + + return BoundingBox(corners); + } + + + + template + std::vector + setup_lexicographic_mask( + const std::array &subcell_indices, + unsigned int dofs_per_line) + { + // This code expands the tensor product space, with `mask` increasing + // in size by a factor `2` for each iteration in `d`. Note how the + // indices in the next line are appended at the end of the vector. + std::vector mask; + mask.reserve(1 << dim); + mask.push_back(0); + unsigned int stride = 1; + for (unsigned int d = 0; d < dim; ++d) + { + const unsigned int previous_mask_size = 1U << d; + AssertDimension(mask.size(), previous_mask_size); + for (unsigned int i = 0; i < previous_mask_size; ++i) + { + mask[i] += subcell_indices[d] * stride; + mask.push_back(stride + mask[i]); + } + stride *= dofs_per_line; + } + return mask; + } } // namespace DiscreteQuadratureGeneratorImplementation } // namespace internal @@ -1627,10 +1795,30 @@ namespace NonMatching + template + void + QuadratureGenerator::clear_quadratures() + { + q_generator.clear_quadratures(); + } + + + template void QuadratureGenerator::generate(const Function & level_set, const BoundingBox &box) + { + clear_quadratures(); + generate_append(level_set, box); + } + + + + template + void + QuadratureGenerator::generate_append(const Function & level_set, + const BoundingBox &box) { Assert(level_set.n_components == 1, ExcMessage( @@ -1638,8 +1826,6 @@ namespace NonMatching " it should have one component.")); Assert(box.volume() > 0, ExcMessage("Incoming box has zero volume.")); - q_generator.clear_quadratures(); - std::vector>> level_sets; level_sets.push_back(level_set); @@ -1703,11 +1889,33 @@ namespace NonMatching + template + void + FaceQuadratureGenerator::clear_quadratures() + { + quadrature_generator.clear_quadratures(); + surface_quadrature.clear(); + } + + + template void FaceQuadratureGenerator::generate(const Function & level_set, const BoundingBox &box, const unsigned int face_index) + { + clear_quadratures(); + generate_append(level_set, box, face_index); + } + + + + template + void + FaceQuadratureGenerator::generate_append(const Function &level_set, + const BoundingBox &box, + const unsigned int face_index) { AssertIndexRange(face_index, GeometryInfo::faces_per_cell); @@ -1727,7 +1935,8 @@ namespace NonMatching // Reuse the lower dimensional QuadratureGenerator on the face. const BoundingBox cross_section = box.cross_section(face_normal_direction); - quadrature_generator.generate(face_restriction, cross_section); + + quadrature_generator.generate_append(face_restriction, cross_section); // We need the dim-dimensional normals of the zero-contour. // Recompute these. @@ -1801,10 +2010,29 @@ namespace NonMatching + void + FaceQuadratureGenerator<1>::clear_quadratures() + { + // do nothing + } + + + void FaceQuadratureGenerator<1>::generate(const Function<1> & level_set, const BoundingBox<1> &box, const unsigned int face_index) + { + clear_quadratures(); + generate_append(level_set, box, face_index); + } + + + + void + FaceQuadratureGenerator<1>::generate_append(const Function<1> & level_set, + const BoundingBox<1> &box, + const unsigned int face_index) { AssertIndexRange(face_index, GeometryInfo<1>::faces_per_cell); @@ -1885,6 +2113,63 @@ namespace NonMatching + template + void + DiscreteQuadratureGenerator::generate_fe_q_iso_q1( + const BoundingBox &unit_box) + { + const unsigned int n_subdivisions = + reference_space_level_set->n_subdivisions(); + + const unsigned int dofs_per_line = n_subdivisions + 1; + + this->clear_quadratures(); + + // The triple nested loop below expands the tensor product of the + // subdivisions, increasing the size of 'all_subcell_indices' by a factor + // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1' + // copies of the data are added at the end of the vector in each + // iteration over d. This code avoids a dimension-dependent + // number of nested loops over each direction, yet having the same + // effect. + std::vector> all_subcell_indices; + all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim)); + all_subcell_indices.emplace_back(); + for (unsigned d = 0; d < dim; ++d) + { + const unsigned int old_size = all_subcell_indices.size(); + for (unsigned int i = 1; i < n_subdivisions; ++i) + { + for (unsigned int j = 0; j < old_size; ++j) + { + std::array next_entry = + all_subcell_indices[j]; + next_entry[d] = i; + all_subcell_indices.push_back(next_entry); + } + } + } + + + for (const auto &subcell_indices : all_subcell_indices) + { + const std::vector lexicographic_mask = + internal::DiscreteQuadratureGeneratorImplementation:: + setup_lexicographic_mask(subcell_indices, dofs_per_line); + + const auto subcell_box = + internal::DiscreteQuadratureGeneratorImplementation:: + create_subcell_box(unit_box, subcell_indices, n_subdivisions); + + reference_space_level_set->set_subcell(lexicographic_mask, subcell_box); + + QuadratureGenerator::generate_append(*reference_space_level_set, + subcell_box); + } + } + + + template void DiscreteQuadratureGenerator::generate( @@ -1896,7 +2181,10 @@ namespace NonMatching reference_space_level_set->set_active_cell(cell); const BoundingBox unit_box = create_unit_bounding_box(); - QuadratureGenerator::generate(*reference_space_level_set, unit_box); + if (reference_space_level_set->is_fe_q_iso_q1()) + generate_fe_q_iso_q1(unit_box); + else + QuadratureGenerator::generate(*reference_space_level_set, unit_box); } @@ -1918,6 +2206,71 @@ namespace NonMatching + template + void + DiscreteFaceQuadratureGenerator::generate_fe_q_iso_q1( + const BoundingBox &unit_box, + unsigned int face_index) + { + const unsigned int n_subdivisions = + reference_space_level_set->n_subdivisions(); + + const unsigned int dofs_per_line = n_subdivisions + 1; + + this->clear_quadratures(); + + const unsigned int coordinate_direction_face = face_index / 2; + const unsigned int coordinate_orientation_face = face_index % 2; + + // The triple nested loop below expands the tensor product of the + // subdivisions, increasing the size of 'all_subcell_indices' by a factor + // 'n_subdivisions' in each iteration. Note how 'n_subdivisions - 1' + // copies of the data are added at the end of the vector in each + // iteration over d. This code avoids a dimension-dependent + // number of nested loops over each direction, yet having the same + // effect. To account for the face all subcell indices in the coordinate of + // the face direction are filled with the appropriate index depending on the + // coordinate orientation. The other two indices (on the face) are + // set accordingly in the nested loops (with a warp around). + std::vector> all_subcell_indices; + all_subcell_indices.reserve(Utilities::pow(n_subdivisions, dim)); + all_subcell_indices.emplace_back(); + all_subcell_indices[0][coordinate_direction_face] = + coordinate_orientation_face == 0 ? 0 : n_subdivisions - 1; + for (unsigned d = 0; d < dim - 1; ++d) + { + const unsigned int old_size = all_subcell_indices.size(); + for (unsigned int i = 1; i < n_subdivisions; ++i) + { + for (unsigned int j = 0; j < old_size; ++j) + { + std::array next_entry = + all_subcell_indices[j]; + next_entry[(coordinate_direction_face + d + 1) % dim] = i; + all_subcell_indices.push_back(next_entry); + } + } + } + + for (const auto &subcell_indices : all_subcell_indices) + { + const std::vector lexicographic_mask = + internal::DiscreteQuadratureGeneratorImplementation:: + setup_lexicographic_mask(subcell_indices, dofs_per_line); + + const auto subcell_box = + internal::DiscreteQuadratureGeneratorImplementation:: + create_subcell_box(unit_box, subcell_indices, n_subdivisions); + + reference_space_level_set->set_subcell(lexicographic_mask, subcell_box); + + FaceQuadratureGenerator::generate_append( + *reference_space_level_set, subcell_box, face_index); + } + } + + + template void DiscreteFaceQuadratureGenerator::generate( @@ -1930,9 +2283,12 @@ namespace NonMatching reference_space_level_set->set_active_cell(cell); const BoundingBox unit_box = create_unit_bounding_box(); - FaceQuadratureGenerator::generate(*reference_space_level_set, - unit_box, - face_index); + if (reference_space_level_set->is_fe_q_iso_q1()) + generate_fe_q_iso_q1(unit_box, face_index); + else + FaceQuadratureGenerator::generate(*reference_space_level_set, + unit_box, + face_index); } } // namespace NonMatching #include "quadrature_generator.inst" diff --git a/tests/non_matching/discrete_quadrature_generator.cc b/tests/non_matching/discrete_quadrature_generator.cc index 87d1be03ce..09b5413372 100644 --- a/tests/non_matching/discrete_quadrature_generator.cc +++ b/tests/non_matching/discrete_quadrature_generator.cc @@ -18,6 +18,7 @@ #include #include +#include #include @@ -35,7 +36,7 @@ template class Test { public: - Test(); + Test(const FE_Poly &fe); void run(); @@ -67,11 +68,15 @@ private: template -Test::Test() +Test::Test(const FE_Poly &fe) : dof_handler(triangulation) { - fe_collection.push_back(FE_Q(1)); - const unsigned int n_quadrature_points = 1; + const FE_Q_iso_Q1 *fe_q_iso_q1 = + dynamic_cast *>(&fe); + const bool is_fe_q_iso_q1 = fe_q_iso_q1 != nullptr; + fe_collection.push_back(fe); + + const unsigned int n_quadrature_points = is_fe_q_iso_q1 ? 1 : fe.get_degree(); q_collection1D.push_back(QGauss<1>(n_quadrature_points)); } @@ -131,8 +136,11 @@ Test::test_discrete_quadrature_generator() q_collection1D, dof_handler, level_set); quadrature_generator.generate(triangulation.begin_active()); + deallog << "inside" << std::endl; print_quadrature(quadrature_generator.get_inside_quadrature()); + deallog << "outside" << std::endl; print_quadrature(quadrature_generator.get_outside_quadrature()); + deallog << "surface" << std::endl; print_surface_quadrature(quadrature_generator.get_surface_quadrature()); } @@ -149,9 +157,12 @@ Test::test_discrete_face_quadrature_generator() for (const auto f : cell->face_indices()) { face_quadrature_generator.generate(cell, f); - + deallog << "face_index = " << f << std::endl; + deallog << "inside" << std::endl; print_quadrature(face_quadrature_generator.get_inside_quadrature()); + deallog << "outside" << std::endl; print_quadrature(face_quadrature_generator.get_outside_quadrature()); + deallog << "surface" << std::endl; print_surface_quadrature( face_quadrature_generator.get_surface_quadrature()); } @@ -161,10 +172,10 @@ Test::test_discrete_face_quadrature_generator() template void -run_test() +run_test(const FE_Poly &fe) { deallog << "dim = " << dim << std::endl; - Test test; + Test test(fe); test.run(); deallog << std::endl; } @@ -176,7 +187,19 @@ main() { initlog(); - run_test<1>(); - run_test<2>(); - run_test<3>(); + std::shared_ptr> fe_1; + std::shared_ptr> fe_2; + std::shared_ptr> fe_3; + fe_1 = std::make_shared>(1); + fe_2 = std::make_shared>(1); + fe_3 = std::make_shared>(1); + run_test<1>(*fe_1); + run_test<2>(*fe_2); + run_test<3>(*fe_3); + fe_1 = std::make_shared>(2); + fe_2 = std::make_shared>(2); + fe_3 = std::make_shared>(2); + run_test<1>(*fe_1); + run_test<2>(*fe_2); + run_test<3>(*fe_3); } diff --git a/tests/non_matching/discrete_quadrature_generator.output b/tests/non_matching/discrete_quadrature_generator.output index 4dfaa1998f..59b2760ca4 100644 --- a/tests/non_matching/discrete_quadrature_generator.output +++ b/tests/non_matching/discrete_quadrature_generator.output @@ -1,40 +1,242 @@ DEAL::dim = 1 +DEAL::inside DEAL::0.375000, 0.750000 +DEAL::outside DEAL::0.875000, 0.250000 +DEAL::surface DEAL::0.750000, 1.00000, 1.00000 +DEAL::face_index = 0 +DEAL::inside DEAL::1.00000 +DEAL::outside +DEAL::surface +DEAL::face_index = 1 +DEAL::inside +DEAL::outside DEAL::1.00000 +DEAL::surface DEAL:: DEAL::dim = 2 +DEAL::inside DEAL::0.375000, 0.500000, 0.750000 +DEAL::outside DEAL::0.875000, 0.500000, 0.250000 +DEAL::surface DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000 +DEAL::face_index = 0 +DEAL::inside DEAL::0.500000, 1.00000 +DEAL::outside +DEAL::surface +DEAL::face_index = 1 +DEAL::inside +DEAL::outside DEAL::0.500000, 1.00000 +DEAL::surface +DEAL::face_index = 2 +DEAL::inside DEAL::0.375000, 0.750000 +DEAL::outside DEAL::0.875000, 0.250000 +DEAL::surface DEAL::0.750000, 1.00000, 1.00000, 0.00000 +DEAL::face_index = 3 +DEAL::inside DEAL::0.375000, 0.750000 +DEAL::outside DEAL::0.875000, 0.250000 +DEAL::surface DEAL::0.750000, 1.00000, 1.00000, 0.00000 DEAL:: DEAL::dim = 3 +DEAL::inside DEAL::0.375000, 0.500000, 0.500000, 0.750000 +DEAL::outside DEAL::0.875000, 0.500000, 0.500000, 0.250000 +DEAL::surface DEAL::0.750000, 0.500000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 0 +DEAL::inside DEAL::0.500000, 0.500000, 1.00000 +DEAL::outside +DEAL::surface +DEAL::face_index = 1 +DEAL::inside +DEAL::outside DEAL::0.500000, 0.500000, 1.00000 +DEAL::surface +DEAL::face_index = 2 +DEAL::inside DEAL::0.500000, 0.375000, 0.750000 +DEAL::outside DEAL::0.500000, 0.875000, 0.250000 +DEAL::surface DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 3 +DEAL::inside DEAL::0.500000, 0.375000, 0.750000 +DEAL::outside DEAL::0.500000, 0.875000, 0.250000 +DEAL::surface DEAL::0.500000, 0.750000, 1.00000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 4 +DEAL::inside DEAL::0.375000, 0.500000, 0.750000 +DEAL::outside DEAL::0.875000, 0.500000, 0.250000 +DEAL::surface DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 5 +DEAL::inside DEAL::0.375000, 0.500000, 0.750000 +DEAL::outside DEAL::0.875000, 0.500000, 0.250000 +DEAL::surface DEAL::0.750000, 0.500000, 1.00000, 1.00000, 0.00000, 0.00000 DEAL:: +DEAL::dim = 1 +DEAL::inside +DEAL::0.250000, 0.500000 +DEAL::0.625000, 0.250000 +DEAL::outside +DEAL::0.875000, 0.250000 +DEAL::surface +DEAL::0.750000, 1.00000, 1.00000 +DEAL::face_index = 0 +DEAL::inside +DEAL::1.00000 +DEAL::outside +DEAL::surface +DEAL::face_index = 1 +DEAL::inside +DEAL::outside +DEAL::1.00000 +DEAL::surface +DEAL:: +DEAL::dim = 2 +DEAL::inside +DEAL::0.250000, 0.250000, 0.250000 +DEAL::0.625000, 0.250000, 0.125000 +DEAL::0.250000, 0.750000, 0.250000 +DEAL::0.625000, 0.750000, 0.125000 +DEAL::outside +DEAL::0.875000, 0.250000, 0.125000 +DEAL::0.875000, 0.750000, 0.125000 +DEAL::surface +DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000 +DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000 +DEAL::face_index = 0 +DEAL::inside +DEAL::0.250000, 0.500000 +DEAL::0.750000, 0.500000 +DEAL::outside +DEAL::surface +DEAL::face_index = 1 +DEAL::inside +DEAL::outside +DEAL::0.250000, 0.500000 +DEAL::0.750000, 0.500000 +DEAL::surface +DEAL::face_index = 2 +DEAL::inside +DEAL::0.250000, 0.500000 +DEAL::0.625000, 0.250000 +DEAL::outside +DEAL::0.875000, 0.250000 +DEAL::surface +DEAL::0.750000, 1.00000, 1.00000, 0.00000 +DEAL::face_index = 3 +DEAL::inside +DEAL::0.250000, 0.500000 +DEAL::0.625000, 0.250000 +DEAL::outside +DEAL::0.875000, 0.250000 +DEAL::surface +DEAL::0.750000, 1.00000, 1.00000, 0.00000 +DEAL:: +DEAL::dim = 3 +DEAL::inside +DEAL::0.250000, 0.250000, 0.250000, 0.125000 +DEAL::0.625000, 0.250000, 0.250000, 0.0625000 +DEAL::0.250000, 0.750000, 0.250000, 0.125000 +DEAL::0.625000, 0.750000, 0.250000, 0.0625000 +DEAL::0.250000, 0.250000, 0.750000, 0.125000 +DEAL::0.625000, 0.250000, 0.750000, 0.0625000 +DEAL::0.250000, 0.750000, 0.750000, 0.125000 +DEAL::0.625000, 0.750000, 0.750000, 0.0625000 +DEAL::outside +DEAL::0.875000, 0.250000, 0.250000, 0.0625000 +DEAL::0.875000, 0.750000, 0.250000, 0.0625000 +DEAL::0.875000, 0.250000, 0.750000, 0.0625000 +DEAL::0.875000, 0.750000, 0.750000, 0.0625000 +DEAL::surface +DEAL::0.750000, 0.250000, 0.250000, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::0.750000, 0.750000, 0.250000, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::0.750000, 0.250000, 0.750000, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::0.750000, 0.750000, 0.750000, 0.250000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 0 +DEAL::inside +DEAL::0.250000, 0.250000, 0.250000 +DEAL::0.750000, 0.250000, 0.250000 +DEAL::0.250000, 0.750000, 0.250000 +DEAL::0.750000, 0.750000, 0.250000 +DEAL::outside +DEAL::surface +DEAL::face_index = 1 +DEAL::inside +DEAL::outside +DEAL::0.250000, 0.250000, 0.250000 +DEAL::0.750000, 0.250000, 0.250000 +DEAL::0.250000, 0.750000, 0.250000 +DEAL::0.750000, 0.750000, 0.250000 +DEAL::surface +DEAL::face_index = 2 +DEAL::inside +DEAL::0.250000, 0.250000, 0.250000 +DEAL::0.750000, 0.250000, 0.250000 +DEAL::0.250000, 0.625000, 0.125000 +DEAL::0.750000, 0.625000, 0.125000 +DEAL::outside +DEAL::0.250000, 0.875000, 0.125000 +DEAL::0.750000, 0.875000, 0.125000 +DEAL::surface +DEAL::0.250000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 3 +DEAL::inside +DEAL::0.250000, 0.250000, 0.250000 +DEAL::0.750000, 0.250000, 0.250000 +DEAL::0.250000, 0.625000, 0.125000 +DEAL::0.750000, 0.625000, 0.125000 +DEAL::outside +DEAL::0.250000, 0.875000, 0.125000 +DEAL::0.750000, 0.875000, 0.125000 +DEAL::surface +DEAL::0.250000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 4 +DEAL::inside +DEAL::0.250000, 0.250000, 0.250000 +DEAL::0.625000, 0.250000, 0.125000 +DEAL::0.250000, 0.750000, 0.250000 +DEAL::0.625000, 0.750000, 0.125000 +DEAL::outside +DEAL::0.875000, 0.250000, 0.125000 +DEAL::0.875000, 0.750000, 0.125000 +DEAL::surface +DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::face_index = 5 +DEAL::inside +DEAL::0.250000, 0.250000, 0.250000 +DEAL::0.625000, 0.250000, 0.125000 +DEAL::0.250000, 0.750000, 0.250000 +DEAL::0.625000, 0.750000, 0.125000 +DEAL::outside +DEAL::0.875000, 0.250000, 0.125000 +DEAL::0.875000, 0.750000, 0.125000 +DEAL::surface +DEAL::0.750000, 0.250000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::0.750000, 0.750000, 0.500000, 1.00000, 0.00000, 0.00000 +DEAL::