From 17025d62db2d8ff13df571d23a0b4b12dc8d3d7c Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 25 Oct 2021 18:29:55 -0600 Subject: [PATCH] Make DataOutFaces work with simplices. --- .../numerics/data_out_dof_data.templates.h | 97 ++++++++++++++++--- tests/fail/step-04-data_out_faces.output | 88 ++++++----------- 2 files changed, 114 insertions(+), 71 deletions(-) diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h index b48c7f7019..ea788c8023 100644 --- a/include/deal.II/numerics/data_out_dof_data.templates.h +++ b/include/deal.II/numerics/data_out_dof_data.templates.h @@ -282,6 +282,7 @@ namespace internal ReferenceCells::Pyramid); })); + // Decide whether we want to work on cell or face FE(Face)Values objects: if (use_face_values == false) { std::unique_ptr> quadrature_simplex; @@ -325,21 +326,23 @@ namespace internal x_fe_values.resize(finite_elements.size()); for (unsigned int i = 0; i < finite_elements.size(); ++i) { - // check if there is a finite element that is equal to the - // present one, then we can re-use the FEValues object + // Check if one of the previous finite elements is equal to the + // present one. If so, re-use the FEValues object. for (unsigned int j = 0; j < i; ++j) if (finite_elements[i].get() == finite_elements[j].get()) { x_fe_values[i] = x_fe_values[j]; break; } + + // If none was found, create an FEValues object: if (x_fe_values[i].get() == nullptr) { dealii::hp::QCollection quadrature; for (unsigned int j = 0; j < finite_elements[i]->size(); ++j) { - const auto reference_cell = + const dealii::ReferenceCell reference_cell = (*finite_elements[i])[j].reference_cell(); if (reference_cell.is_hyper_cube()) @@ -379,31 +382,84 @@ namespace internal } else // build FEFaceValues objects instead { - dealii::hp::QCollection quadrature( - QIterated(QTrapezoid<1>(), n_subdivisions)); + std::unique_ptr> quadrature_simplex; + std::unique_ptr> quadrature_hypercube; + + // See whether we need simplex or hypercube quadrature formulas. + // This is only an issue in 3d (in 2d every face integral is just + // a line integral, so we can deal with that via the usual hypercube + // quadrature rule). + if ((dim == 3) && + (needs_simplex_setup || needs_pyramid_setup || needs_wedge_setup)) + { + quadrature_simplex = std::make_unique>( + generate_simplex_evaluation_points(n_subdivisions)); + } + + if ((dim < 3) || (needs_hypercube_setup || needs_pyramid_setup || + needs_wedge_setup)) + { + quadrature_hypercube = + std::make_unique>(QTrapezoid<1>(), + n_subdivisions); + } x_fe_face_values.resize(finite_elements.size()); for (unsigned int i = 0; i < finite_elements.size(); ++i) { - // check if there is a finite element that is equal to the - // present one, then we can re-use the FEValues object + // Check if one of the previous finite elements is equal to the + // present one. If so, re-use the FEValues object. for (unsigned int j = 0; j < i; ++j) if (finite_elements[i].get() == finite_elements[j].get()) { x_fe_face_values[i] = x_fe_face_values[j]; break; } + + // If none was found, create an FEFaceValues object: if (x_fe_face_values[i].get() == nullptr) - x_fe_face_values[i] = - std::make_shared>( - mapping_collection, - *finite_elements[i], - quadrature, - update_flags); + { + dealii::hp::QCollection quadrature; + + for (unsigned int j = 0; j < finite_elements[i]->size(); ++j) + { + const dealii::ReferenceCell reference_cell = + (*finite_elements[i])[j].reference_cell(); + + // In 1d/2d and for hypercube/wedge/pyramid elements, we + // need hypercube quadratures. + if ((dim < 3) || + (reference_cell.is_hyper_cube() || + (reference_cell == dealii::ReferenceCells::Wedge) || + (reference_cell == dealii::ReferenceCells::Pyramid))) + quadrature.push_back(*quadrature_hypercube); + + // In 3d, if the element is for simplex/wedge/pyramid + // cells, then we also need simplex quadratures + if ((dim == 3) && + (reference_cell.is_simplex() || + (reference_cell == dealii::ReferenceCells::Wedge) || + (reference_cell == dealii::ReferenceCells::Pyramid))) + quadrature.push_back(*quadrature_simplex); + } + + x_fe_face_values[i] = + std::make_shared>( + mapping_collection, + *finite_elements[i], + quadrature, + update_flags); + } } // Return maximal number of evaluation points: - return quadrature[0].size(); + return std::max( + {(dim == 3) && (needs_simplex_setup || needs_pyramid_setup || + needs_wedge_setup) ? + quadrature_simplex->size() : + 0, + (dim < 3) || needs_hypercube_setup ? quadrature_hypercube->size() : + 0}); } } @@ -496,11 +552,14 @@ namespace internal { if (cell->is_active()) { - typename DoFHandler::active_cell_iterator + const typename DoFHandler::active_cell_iterator dh_cell(&cell->get_triangulation(), cell->level(), cell->index(), dof_data[dataset]->dof_handler); + + // Check whether we need cell or face FEValues objects by + // testing which of the two arrays actually has any content. if (x_fe_values.empty()) { AssertIndexRange(face, GeometryInfo::faces_per_cell); @@ -513,6 +572,11 @@ namespace internal x_fe_values[dataset]->reinit(cell); } } + + // If there is are no DoF-associated data (just cell-associated ones), + // then the loop above will not execute any iterations. In that case, + // do the initialization for the first FE(Face)Values object by + // hand, using only the (triangulation) cell without a DoFHandler. if (dof_data.empty()) { if (x_fe_values.empty()) @@ -533,6 +597,9 @@ namespace internal const unsigned int dataset) const { AssertIndexRange(dataset, finite_elements.size()); + + // Check whether we need cell or face FEValues objects by testing + // which of the two arrays actually has any content. if (x_fe_values.empty()) return x_fe_face_values[dataset]->get_present_fe_values(); else diff --git a/tests/fail/step-04-data_out_faces.output b/tests/fail/step-04-data_out_faces.output index 484c320883..a795ce8626 100644 --- a/tests/fail/step-04-data_out_faces.output +++ b/tests/fail/step-04-data_out_faces.output @@ -1,4 +1,4 @@ - +JobId strange.math.colostate.edu Tue Oct 26 17:19:44 2021 DEAL::Solving problem in 2 space dimensions. DEAL:: Number of active cells: 8 DEAL:: Total number of cells: 8 @@ -7,7 +7,7 @@ DEAL:cg::Starting value 4.18267 DEAL:cg::Convergence step 1 value 0.00000 DEAL:: 1 CG iterations needed to obtain convergence. # vtk DataFile Version 3.0 -#This file was generated +#This file was generated by the deal.II library on 2021/10/26 at 17:19:44 ASCII DATASET UNSTRUCTURED_GRID @@ -53,137 +53,113 @@ DEAL:cg::Starting value 0.00000 DEAL:cg::Convergence step 0 value 0.00000 DEAL:: 0 CG iterations needed to obtain convergence. # vtk DataFile Version 3.0 -#This file was generated +#This file was generated by the deal.II library on 2021/10/26 at 17:19:44 ASCII DATASET UNSTRUCTURED_GRID -POINTS 96 double +POINTS 72 double -1.00000 -1.00000 -1.00000 1.00000 -1.00000 -1.00000 0.00000 0.00000 -1.00000 -2.00000 0.00000 -1.00000 1.00000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 0.00000 -1.00000 0.00000 --2.00000 -1.00000 0.00000 -1.00000 1.00000 -1.00000 1.00000 1.00000 -1.00000 0.00000 1.00000 0.00000 -2.00000 1.00000 0.00000 1.00000 1.00000 -1.00000 -1.00000 1.00000 -1.00000 0.00000 0.00000 -1.00000 --2.00000 0.00000 -1.00000 1.00000 1.00000 1.00000 -1.00000 1.00000 1.00000 0.00000 1.00000 0.00000 --2.00000 1.00000 0.00000 -1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.00000 0.00000 1.00000 -2.00000 0.00000 1.00000 1.00000 -1.00000 1.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 --2.00000 0.00000 1.00000 -1.00000 -1.00000 1.00000 1.00000 -1.00000 1.00000 0.00000 -1.00000 0.00000 -2.00000 -1.00000 0.00000 -1.00000 -1.00000 -1.00000 -1.00000 1.00000 -1.00000 -1.00000 0.00000 0.00000 --1.00000 2.00000 0.00000 -1.00000 1.00000 -1.00000 -1.00000 -1.00000 -1.00000 0.00000 0.00000 -1.00000 -0.00000 -2.00000 -1.00000 -1.00000 -1.00000 1.00000 -1.00000 1.00000 1.00000 0.00000 0.00000 1.00000 -0.00000 2.00000 1.00000 -1.00000 1.00000 1.00000 -1.00000 -1.00000 1.00000 -1.00000 0.00000 0.00000 --1.00000 -2.00000 0.00000 1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 1.00000 1.00000 1.00000 -0.00000 2.00000 1.00000 1.00000 -1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.00000 0.00000 -1.00000 2.00000 0.00000 1.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 1.00000 -1.00000 -1.00000 2.00000 0.00000 1.00000 -1.00000 -1.00000 1.00000 1.00000 -1.00000 0.00000 0.00000 -1.00000 -0.00000 2.00000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 0.00000 0.00000 -1.00000 -1.00000 1.00000 --1.00000 0.00000 2.00000 -1.00000 -1.00000 -1.00000 -1.00000 -1.00000 1.00000 0.00000 -1.00000 0.00000 -0.00000 -1.00000 2.00000 1.00000 -1.00000 -1.00000 1.00000 -1.00000 1.00000 1.00000 0.00000 0.00000 -1.00000 0.00000 2.00000 1.00000 -1.00000 1.00000 1.00000 -1.00000 -1.00000 0.00000 -1.00000 0.00000 -0.00000 -1.00000 -2.00000 1.00000 1.00000 -1.00000 1.00000 1.00000 1.00000 0.00000 1.00000 0.00000 -0.00000 1.00000 2.00000 1.00000 1.00000 1.00000 1.00000 1.00000 -1.00000 1.00000 0.00000 0.00000 -1.00000 0.00000 -2.00000 -1.00000 1.00000 -1.00000 -1.00000 1.00000 1.00000 -1.00000 0.00000 0.00000 --1.00000 0.00000 2.00000 -1.00000 1.00000 1.00000 -1.00000 1.00000 -1.00000 0.00000 1.00000 0.00000 -0.00000 1.00000 -2.00000 -CELLS 24 120 -4 0 1 3 2 -4 4 5 7 6 -4 8 9 11 10 -4 12 13 15 14 -4 16 17 19 18 -4 20 21 23 22 -4 24 25 27 26 -4 28 29 31 30 -4 32 33 35 34 -4 36 37 39 38 -4 40 41 43 42 -4 44 45 47 46 -4 48 49 51 50 -4 52 53 55 54 -4 56 57 59 58 -4 60 61 63 62 -4 64 65 67 66 -4 68 69 71 70 -4 72 73 75 74 -4 76 77 79 78 -4 80 81 83 82 -4 84 85 87 86 -4 88 89 91 90 -4 92 93 95 94 +CELLS 24 96 + 3 0 1 2 + 3 3 4 5 + 3 6 7 8 + 3 9 10 11 + 3 12 13 14 + 3 15 16 17 + 3 18 19 20 + 3 21 22 23 + 3 24 25 26 + 3 27 28 29 + 3 30 31 32 + 3 33 34 35 + 3 36 37 38 + 3 39 40 41 + 3 42 43 44 + 3 45 46 47 + 3 48 49 50 + 3 51 52 53 + 3 54 55 56 + 3 57 58 59 + 3 60 61 62 + 3 63 64 65 + 3 66 67 68 + 3 69 70 71 CELL_TYPES 24 - 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 -POINT_DATA 96 + 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 +POINT_DATA 72 SCALARS solution double 1 LOOKUP_TABLE default -3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 +3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 -- 2.39.5