From: David Wells Date: Fri, 30 Jun 2023 22:21:26 +0000 (-0400) Subject: Add FE_SimplexP(3) and FE_SimplexDGP(3). X-Git-Tag: relicensing~395^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d792d234b3d71a2f957a8d374e00fcfd0327f016;p=dealii.git Add FE_SimplexP(3) and FE_SimplexDGP(3). --- diff --git a/doc/news/changes/major/20231011DavidWellsPeterMunch b/doc/news/changes/major/20231011DavidWellsPeterMunch new file mode 100644 index 0000000000..58709a9890 --- /dev/null +++ b/doc/news/changes/major/20231011DavidWellsPeterMunch @@ -0,0 +1,3 @@ +New: FE_SimplexP and FE_SimplexDGP now support cubic elements. +
+(David Wells and Peter Munch, 2023/10/11) diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc index 32f42a2167..767d5b1ab5 100644 --- a/source/fe/fe_simplex_p.cc +++ b/source/fe/fe_simplex_p.cc @@ -28,6 +28,31 @@ DEAL_II_NAMESPACE_OPEN namespace { + // TODO: replace this with QProjector once QProjector learns how to project + // quadrature points onto separate faces of simplices + template + Point + face_midpoint(const unsigned int face_no) + { + const auto reference_cell = ReferenceCells::get_simplex(); + const auto face_reference_cell = + reference_cell.face_reference_cell(face_no); + + Point midpoint; + for (const auto face_vertex_no : face_reference_cell.vertex_indices()) + { + const auto vertex_no = reference_cell.face_to_cell_vertices( + face_no, + face_vertex_no, + ReferenceCell::default_combined_face_orientation()); + + midpoint += reference_cell.template vertex(vertex_no); + } + + midpoint /= reference_cell.face_reference_cell(0).n_vertices(); + return midpoint; + } + /** * Helper function to set up the dpo vector of FE_SimplexP for a given @p dim and * @p degree. @@ -44,6 +69,8 @@ namespace return {1, 0}; case 2: return {1, 1}; + case 3: + return {1, 2}; default: Assert(false, ExcNotImplemented()); } @@ -54,6 +81,8 @@ namespace return {1, 0, 0}; case 2: return {1, 1, 0}; + case 3: + return {1, 2, 1}; default: Assert(false, ExcNotImplemented()); } @@ -64,6 +93,8 @@ namespace return {1, 0, 0, 0}; case 2: return {1, 1, 0, 0}; + case 3: + return {1, 2, 1, 0}; default: Assert(false, ExcNotImplemented()); } @@ -83,15 +114,11 @@ namespace std::vector> unit_support_points_fe_p(const unsigned int degree) { + Assert(dim != 0, ExcInternalError()); + Assert(degree <= 3, ExcNotImplemented()); std::vector> unit_points; - // If we do dim - 1 we can get here in dim = 0: - if (dim == 0) - { - unit_points.emplace_back(); - return unit_points; - } + const auto reference_cell = ReferenceCells::get_simplex(); - const auto reference_cell = ReferenceCells::get_simplex(); // Piecewise constants are a special case: use a support point at the // centroid and only the centroid if (degree == 0) @@ -100,23 +127,45 @@ namespace return unit_points; } - Assert(degree <= 2, ExcNotImplemented()); - for (const auto &vertex_no : reference_cell.vertex_indices()) - unit_points.emplace_back(reference_cell.template vertex(vertex_no)); - - if (degree == 2) - for (const auto &line_no : reference_cell.line_indices()) - { - const auto v0 = reference_cell.template vertex( - reference_cell.line_to_cell_vertices(line_no, 0)); - const auto v1 = reference_cell.template vertex( - reference_cell.line_to_cell_vertices(line_no, 1)); - unit_points.emplace_back((v0 + v1) / 2.0); - } + // otherwise write everything as linear combinations of vertices + const auto dpo = get_dpo_vector_fe_p(dim, degree); + Assert(dpo.size() == dim + 1, ExcInternalError()); + Assert(dpo[0] == 1, ExcNotImplemented()); + // vertices: + for (const unsigned int d : reference_cell.vertex_indices()) + unit_points.push_back(reference_cell.template vertex(d)); + // lines: + for (const unsigned int l : reference_cell.line_indices()) + { + const Point p0 = + unit_points[reference_cell.line_to_cell_vertices(l, 0)]; + const Point p1 = + unit_points[reference_cell.line_to_cell_vertices(l, 1)]; + for (unsigned int p = 0; p < dpo[1]; ++p) + unit_points.push_back((double(dpo[1] - p) / (dpo[1] + 1)) * p0 + + (double(p + 1) / (dpo[1] + 1)) * p1); + } + // quads: + if (dim > 1 && dpo[2] > 0) + { + Assert(dpo[2] == 1, ExcNotImplemented()); + if (dim == 2) + unit_points.push_back(reference_cell.template barycenter()); + if (dim == 3) + for (const unsigned int f : reference_cell.face_indices()) + unit_points.push_back(face_midpoint(f)); + } return unit_points; } + template <> + std::vector> + unit_support_points_fe_p(const unsigned int /*degree*/) + { + return {Point<0>()}; + } + /** * Set up a vector that contains the unit (reference) cell's faces support * points for FE_SimplexP and sufficiently similar elements. @@ -166,9 +215,8 @@ namespace FullMatrix constraints_fe_p<2>(const unsigned int degree) { - const unsigned int dim = 2; - - Assert(degree <= 2, ExcNotImplemented()); + constexpr int dim = 2; + Assert(degree <= 3, ExcNotImplemented()); // the following implements the 2d case // (the 3d case is not implemented yet) @@ -575,7 +623,14 @@ FE_SimplexP::FE_SimplexP(const unsigned int degree) unit_support_points_fe_p(degree), unit_face_support_points_fe_p(degree, FiniteElementData::H1), constraints_fe_p(degree)) -{} +{ + if (degree > 2) + for (unsigned int i = 0; i < this->n_dofs_per_line(); ++i) + this->adjust_line_dof_index_for_line_orientation_table[i] = + this->n_dofs_per_line() - 1 - i - i; + // We do not support multiple DoFs per quad yet + Assert(degree <= 3, ExcNotImplemented()); +} diff --git a/tests/simplex/simplex_project_cg.cc b/tests/simplex/simplex_project_cg.cc index 31132a4752..b0377afbf6 100644 --- a/tests/simplex/simplex_project_cg.cc +++ b/tests/simplex/simplex_project_cg.cc @@ -84,12 +84,13 @@ test(const unsigned int degree) { #ifdef SIMPLEX FE_SimplexP fe(degree); - QGaussSimplex quadrature(degree + 2); + QGaussSimplex quadrature(degree + 1); #else FE_Q fe(degree); - QGauss quadrature(degree + 2); + QGauss quadrature(degree + 1); #endif deallog << "FE = " << fe.get_name() << std::endl; + deallog << std::setprecision(4); double previous_error = 1.0; @@ -116,65 +117,52 @@ test(const unsigned int degree) reference_cell.template get_default_linear_mapping(); dummy.close(); - const bool l2_projection = false; + SparsityPattern sparsity_pattern; + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + } - if (l2_projection) - { - VectorTools::project( - mapping, dof_handler, dummy, quadrature, function, solution); - } - else - { - SparsityPattern sparsity_pattern; - { - DynamicSparsityPattern dsp(dof_handler.n_dofs(), - dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dsp); - sparsity_pattern.copy_from(dsp); - } - - SparseMatrix h1_matrix(sparsity_pattern); - SparseMatrix laplace_matrix(sparsity_pattern); - - MatrixCreator::create_mass_matrix(mapping, - dof_handler, - quadrature, - h1_matrix); - MatrixCreator::create_laplace_matrix(mapping, - dof_handler, - quadrature, - laplace_matrix); - - h1_matrix.add(0.0, laplace_matrix); - - Vector rhs(solution.size()); - VectorTools::create_right_hand_side( - mapping, dof_handler, quadrature, ForcingH1(), rhs); - - SolverControl solver_control(1000, - 1e-12 * rhs.l2_norm(), - false, - false); - SolverCG> cg(solver_control); - - cg.solve(h1_matrix, solution, rhs, PreconditionIdentity()); - } + SparseMatrix h1_matrix(sparsity_pattern); + SparseMatrix laplace_matrix(sparsity_pattern); + + MatrixCreator::create_mass_matrix(mapping, + dof_handler, + quadrature, + h1_matrix); + MatrixCreator::create_laplace_matrix(mapping, + dof_handler, + quadrature, + laplace_matrix); + + h1_matrix.add(0.0, laplace_matrix); + + Vector rhs(solution.size()); + VectorTools::create_right_hand_side( + mapping, dof_handler, quadrature, ForcingH1(), rhs); + + SolverControl solver_control(1000, 1e-12 * rhs.l2_norm(), false, false); + SolverCG> cg(solver_control); + + cg.solve(h1_matrix, solution, rhs, PreconditionIdentity()); VectorTools::integrate_difference(mapping, dof_handler, solution, function, cell_errors, - Quadrature( - fe.get_unit_support_points()), - VectorTools::Linfty_norm); - - const double max_error = - *std::max_element(cell_errors.begin(), cell_errors.end()); - deallog << "max error = " << max_error << std::endl; - if (max_error != 0.0) - deallog << "ratio = " << previous_error / max_error << std::endl; - previous_error = max_error; + quadrature, + VectorTools::L2_norm); + const double L2_error = + VectorTools::compute_global_error(tria, + cell_errors, + VectorTools::L2_norm); + + deallog << "L2 error = " << L2_error << std::endl; + if (L2_error != 0.0) + deallog << "ratio = " << previous_error / L2_error << std::endl; + previous_error = L2_error; #if 0 if (dim == 2) @@ -199,7 +187,9 @@ main() test<2>(1); test<2>(2); + test<2>(3); test<3>(1); test<3>(2); + test<3>(3); } diff --git a/tests/simplex/simplex_project_cg.output b/tests/simplex/simplex_project_cg.output index 287c3090a0..dfd98e33ac 100644 --- a/tests/simplex/simplex_project_cg.output +++ b/tests/simplex/simplex_project_cg.output @@ -1,37 +1,55 @@ DEAL::FE = FE_SimplexP<2>(1) -DEAL::max error = 0.332362 -DEAL::ratio = 3.00877 -DEAL::max error = 0.100426 -DEAL::ratio = 3.30953 -DEAL::max error = 0.0271946 -DEAL::ratio = 3.69285 -DEAL::max error = 0.00747939 -DEAL::ratio = 3.63594 +DEAL::L2 error = 0.05224 +DEAL::ratio = 19.14 +DEAL::L2 error = 0.02949 +DEAL::ratio = 1.772 +DEAL::L2 error = 0.007768 +DEAL::ratio = 3.796 +DEAL::L2 error = 0.001951 +DEAL::ratio = 3.981 DEAL::FE = FE_SimplexP<2>(2) -DEAL::max error = 0.0430400 -DEAL::ratio = 23.2342 -DEAL::max error = 0.0116055 -DEAL::ratio = 3.70859 -DEAL::max error = 0.00232162 -DEAL::ratio = 4.99887 -DEAL::max error = 0.000336599 -DEAL::ratio = 6.89729 +DEAL::L2 error = 0.004078 +DEAL::ratio = 245.2 +DEAL::L2 error = 0.002064 +DEAL::ratio = 1.975 +DEAL::L2 error = 0.0003512 +DEAL::ratio = 5.879 +DEAL::L2 error = 5.040e-05 +DEAL::ratio = 6.967 +DEAL::FE = FE_SimplexP<2>(3) +DEAL::L2 error = 0.004061 +DEAL::ratio = 246.3 +DEAL::L2 error = 0.0002464 +DEAL::ratio = 16.48 +DEAL::L2 error = 1.603e-05 +DEAL::ratio = 15.38 +DEAL::L2 error = 1.009e-06 +DEAL::ratio = 15.88 DEAL::FE = FE_SimplexP<3>(1) -DEAL::max error = 0.620494 -DEAL::ratio = 1.61162 -DEAL::max error = 0.259138 -DEAL::ratio = 2.39445 -DEAL::max error = 0.0619310 -DEAL::ratio = 4.18431 -DEAL::max error = 0.0151022 -DEAL::ratio = 4.10080 +DEAL::L2 error = 0.09826 +DEAL::ratio = 10.18 +DEAL::L2 error = 0.02938 +DEAL::ratio = 3.344 +DEAL::L2 error = 0.007387 +DEAL::ratio = 3.977 +DEAL::L2 error = 0.001763 +DEAL::ratio = 4.191 DEAL::FE = FE_SimplexP<3>(2) -DEAL::max error = 0.0550999 -DEAL::ratio = 18.1489 -DEAL::max error = 0.0466852 -DEAL::ratio = 1.18024 -DEAL::max error = 0.00502563 -DEAL::ratio = 9.28942 -DEAL::max error = 0.000810078 -DEAL::ratio = 6.20388 +DEAL::L2 error = 0.01087 +DEAL::ratio = 91.95 +DEAL::L2 error = 0.004092 +DEAL::ratio = 2.657 +DEAL::L2 error = 0.0007005 +DEAL::ratio = 5.842 +DEAL::L2 error = 0.0001018 +DEAL::ratio = 6.880 +DEAL::FE = FE_SimplexP<3>(3) +DEAL::L2 error = 0.006416 +DEAL::ratio = 155.9 +DEAL::L2 error = 0.0003557 +DEAL::ratio = 18.03 +DEAL::L2 error = 2.324e-05 +DEAL::ratio = 15.31 +DEAL::L2 error = 1.447e-06 +DEAL::ratio = 16.06 diff --git a/tests/simplex/simplex_project_dg.cc b/tests/simplex/simplex_project_dg.cc index 5edb6fae80..4d277490c6 100644 --- a/tests/simplex/simplex_project_dg.cc +++ b/tests/simplex/simplex_project_dg.cc @@ -20,18 +20,29 @@ #include #include +#include +#include #include +#include #include #include #include +#include +#include +#include +#include +#include +#include #include #include +#include #include #include +#include #include "../tests.h" @@ -41,6 +52,7 @@ test(const unsigned int degree) { FE_SimplexDGP fe(degree); deallog << "FE = " << fe.get_name() << std::endl; + deallog << std::setprecision(4); QGaussSimplex quadrature(degree + 1); double previous_error = 1.0; @@ -63,24 +75,44 @@ test(const unsigned int degree) const auto &mapping = reference_cell.template get_default_linear_mapping(); dummy.close(); - VectorTools::project( - mapping, dof_handler, dummy, quadrature, function, solution); + + SparsityPattern sparsity_pattern; + { + DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + } + + SparseMatrix mass_matrix(sparsity_pattern); + MatrixCreator::create_mass_matrix(mapping, + dof_handler, + quadrature, + mass_matrix); + Vector rhs(solution.size()); + VectorTools::create_right_hand_side( + mapping, dof_handler, quadrature, function, rhs); + + SolverControl solver_control(1000, 1e-12 * rhs.l2_norm(), false, false); + SolverCG> cg(solver_control); + + cg.solve(mass_matrix, solution, rhs, PreconditionIdentity()); VectorTools::integrate_difference(mapping, dof_handler, solution, function, cell_errors, - Quadrature( - fe.get_unit_support_points()), - VectorTools::Linfty_norm); - - const double max_error = - *std::max_element(cell_errors.begin(), cell_errors.end()); - deallog << "max error = " << max_error << std::endl; - if (max_error != 0.0) - deallog << "ratio = " << previous_error / max_error << std::endl; - previous_error = max_error; + quadrature, + VectorTools::L2_norm); + const double L2_error = + VectorTools::compute_global_error(tria, + cell_errors, + VectorTools::L2_norm); + + deallog << "L2 error = " << L2_error << std::endl; + if (L2_error != 0.0) + deallog << "ratio = " << previous_error / L2_error << std::endl; + previous_error = L2_error; #if 0 if (dim == 2) @@ -105,7 +137,9 @@ main() test<2>(1); test<2>(2); + test<2>(3); test<3>(1); test<3>(2); + test<3>(3); } diff --git a/tests/simplex/simplex_project_dg.output b/tests/simplex/simplex_project_dg.output index 4bb911968d..1b38184f1a 100644 --- a/tests/simplex/simplex_project_dg.output +++ b/tests/simplex/simplex_project_dg.output @@ -1,37 +1,55 @@ DEAL::FE = FE_SimplexDGP<2>(1) -DEAL::max error = 0.113889 -DEAL::ratio = 8.78048 -DEAL::max error = 0.0302392 -DEAL::ratio = 3.76627 -DEAL::max error = 0.00767275 -DEAL::ratio = 3.94112 -DEAL::max error = 0.00192529 -DEAL::ratio = 3.98525 +DEAL::L2 error = 0.009794 +DEAL::ratio = 102.1 +DEAL::L2 error = 0.002479 +DEAL::ratio = 3.952 +DEAL::L2 error = 0.0006215 +DEAL::ratio = 3.988 +DEAL::L2 error = 0.0001555 +DEAL::ratio = 3.997 DEAL::FE = FE_SimplexDGP<2>(2) -DEAL::max error = 0.0174164 -DEAL::ratio = 57.4171 -DEAL::max error = 0.00227404 -DEAL::ratio = 7.65880 -DEAL::max error = 0.000287343 -DEAL::ratio = 7.91404 -DEAL::max error = 3.60148e-05 -DEAL::ratio = 7.97847 +DEAL::L2 error = 0.0002017 +DEAL::ratio = 4957. +DEAL::L2 error = 2.334e-05 +DEAL::ratio = 8.644 +DEAL::L2 error = 2.854e-06 +DEAL::ratio = 8.178 +DEAL::L2 error = 3.547e-07 +DEAL::ratio = 8.046 +DEAL::FE = FE_SimplexDGP<2>(3) +DEAL::L2 error = 0.0001932 +DEAL::ratio = 5175. +DEAL::L2 error = 1.224e-05 +DEAL::ratio = 15.79 +DEAL::L2 error = 7.677e-07 +DEAL::ratio = 15.95 +DEAL::L2 error = 4.802e-08 +DEAL::ratio = 15.99 DEAL::FE = FE_SimplexDGP<3>(1) -DEAL::max error = 0.208217 -DEAL::ratio = 4.80269 -DEAL::max error = 0.0629409 -DEAL::ratio = 3.30813 -DEAL::max error = 0.0164605 -DEAL::ratio = 3.82376 -DEAL::max error = 0.00416117 -DEAL::ratio = 3.95573 +DEAL::L2 error = 0.01241 +DEAL::ratio = 80.59 +DEAL::L2 error = 0.003150 +DEAL::ratio = 3.939 +DEAL::L2 error = 0.0007879 +DEAL::ratio = 3.998 +DEAL::L2 error = 0.0001968 +DEAL::ratio = 4.004 DEAL::FE = FE_SimplexDGP<3>(2) -DEAL::max error = 0.0415553 -DEAL::ratio = 24.0643 -DEAL::max error = 0.00641441 -DEAL::ratio = 6.47844 -DEAL::max error = 0.000842371 -DEAL::ratio = 7.61471 -DEAL::max error = 0.000106584 -DEAL::ratio = 7.90336 +DEAL::L2 error = 0.001633 +DEAL::ratio = 612.3 +DEAL::L2 error = 0.0002072 +DEAL::ratio = 7.883 +DEAL::L2 error = 2.599e-05 +DEAL::ratio = 7.971 +DEAL::L2 error = 3.252e-06 +DEAL::ratio = 7.993 +DEAL::FE = FE_SimplexDGP<3>(3) +DEAL::L2 error = 0.0002064 +DEAL::ratio = 4845. +DEAL::L2 error = 1.309e-05 +DEAL::ratio = 15.77 +DEAL::L2 error = 8.214e-07 +DEAL::ratio = 15.94 +DEAL::L2 error = 5.138e-08 +DEAL::ratio = 15.99