From c67d48d2239f1dcbe410544efd9107a8d54ffb7c Mon Sep 17 00:00:00 2001 From: David Wells Date: Tue, 15 Oct 2024 14:25:13 -0400 Subject: [PATCH] FE_SimplexP: support hanging nodes with cubic elements. This more-or-less just copies the same code used to place interpolation points on child lines we use in FE_Q_Base. --- doc/news/changes/minor/20241015DavidWells | 3 ++ source/fe/fe_simplex_p.cc | 18 ++++++----- tests/simplex/simplex_project_cg.cc | 39 +++++++++++++++-------- tests/simplex/simplex_project_cg.output | 24 ++++++++++++++ 4 files changed, 62 insertions(+), 22 deletions(-) create mode 100644 doc/news/changes/minor/20241015DavidWells diff --git a/doc/news/changes/minor/20241015DavidWells b/doc/news/changes/minor/20241015DavidWells new file mode 100644 index 0000000000..2eca18177f --- /dev/null +++ b/doc/news/changes/minor/20241015DavidWells @@ -0,0 +1,3 @@ +Fixed: FESystem now works correctly with cubic FE_SimplexP elements. +
+(David Wells, 2024/10/15) diff --git a/source/fe/fe_simplex_p.cc b/source/fe/fe_simplex_p.cc index d3cede7d32..1522a51420 100644 --- a/source/fe/fe_simplex_p.cc +++ b/source/fe/fe_simplex_p.cc @@ -218,7 +218,6 @@ namespace constraints_fe_p<2>(const unsigned int degree) { constexpr int dim = 2; - Assert(degree <= 3, ExcNotImplemented()); // the following implements the 2d case // (the 3d case is not implemented yet) @@ -229,13 +228,16 @@ namespace std::vector> constraint_points; // midpoint constraint_points.emplace_back(0.5); - if (degree == 2) - { - // midpoint on subface 0 - constraint_points.emplace_back(0.25); - // midpoint on subface 1 - constraint_points.emplace_back(0.75); - } + // subface 0 + for (unsigned int i = 1; i < degree; ++i) + constraint_points.push_back( + GeometryInfo::child_to_cell_coordinates( + Point(i / double(degree)), 0)); + // subface 1 + for (unsigned int i = 1; i < degree; ++i) + constraint_points.push_back( + GeometryInfo::child_to_cell_coordinates( + Point(i / double(degree)), 1)); // Now construct relation between destination (child) and source (mother) // dofs. diff --git a/tests/simplex/simplex_project_cg.cc b/tests/simplex/simplex_project_cg.cc index 710e2f610a..f79fd1d25a 100644 --- a/tests/simplex/simplex_project_cg.cc +++ b/tests/simplex/simplex_project_cg.cc @@ -79,7 +79,7 @@ public: template void -test(const unsigned int degree) +test(const unsigned int degree, const bool use_amr = false) { #ifdef SIMPLEX FE_SimplexP fe(degree); @@ -90,10 +90,12 @@ test(const unsigned int degree) #endif deallog << "FE = " << fe.get_name() << std::endl; deallog << std::setprecision(4); + if (use_amr) + deallog << "using AMR" << std::endl; double previous_error = 1.0; - for (unsigned int r = 0; r < 4; ++r) + for (unsigned int r = 0; r < (use_amr ? 3 : 4); ++r) { Triangulation tria_hex, tria; GridGenerator::hyper_cube(tria_hex); @@ -103,6 +105,13 @@ test(const unsigned int degree) #else tria.copy_triangulation(tria_hex); #endif + if (use_amr) + { + for (auto &cell : tria.active_cell_iterators()) + if (cell->index() % 3 == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } ReferenceCell reference_cell = tria.begin_active()->reference_cell(); DoFHandler dof_handler(tria); @@ -111,40 +120,38 @@ test(const unsigned int degree) Vector cell_errors(tria.n_active_cells()); Vector solution(dof_handler.n_dofs()); Solution function; - AffineConstraints dummy; + AffineConstraints constraints; const auto &mapping = reference_cell.template get_default_linear_mapping(); - dummy.close(); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); SparsityPattern sparsity_pattern; { DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs()); - DoFTools::make_sparsity_pattern(dof_handler, dsp); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, true); 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); + MatrixCreator::create_mass_matrix( + mapping, dof_handler, quadrature, h1_matrix, {}, constraints); + MatrixCreator::create_laplace_matrix( + mapping, dof_handler, quadrature, laplace_matrix, {}, constraints); h1_matrix.add(0.0, laplace_matrix); Vector rhs(solution.size()); VectorTools::create_right_hand_side( - mapping, dof_handler, quadrature, ForcingH1(), rhs); + mapping, dof_handler, quadrature, ForcingH1(), rhs, constraints); SolverControl solver_control(1000, 1e-12 * rhs.l2_norm(), false, false); SolverCG> cg(solver_control); cg.solve(h1_matrix, solution, rhs, PreconditionIdentity()); + constraints.distribute(solution); VectorTools::integrate_difference(mapping, dof_handler, @@ -191,4 +198,8 @@ main() test<3>(1); test<3>(2); test<3>(3); + + test<2>(1, true); + test<2>(2, true); + test<2>(3, true); } diff --git a/tests/simplex/simplex_project_cg.output b/tests/simplex/simplex_project_cg.output index dfd98e33ac..97c1f6475b 100644 --- a/tests/simplex/simplex_project_cg.output +++ b/tests/simplex/simplex_project_cg.output @@ -53,3 +53,27 @@ DEAL::L2 error = 2.324e-05 DEAL::ratio = 15.31 DEAL::L2 error = 1.447e-06 DEAL::ratio = 16.06 +DEAL::FE = FE_SimplexP<2>(1) +DEAL::using AMR +DEAL::L2 error = 0.07022 +DEAL::ratio = 14.24 +DEAL::L2 error = 0.03106 +DEAL::ratio = 2.261 +DEAL::L2 error = 0.008068 +DEAL::ratio = 3.850 +DEAL::FE = FE_SimplexP<2>(2) +DEAL::using AMR +DEAL::L2 error = 0.004448 +DEAL::ratio = 224.8 +DEAL::L2 error = 0.002054 +DEAL::ratio = 2.165 +DEAL::L2 error = 0.0003406 +DEAL::ratio = 6.030 +DEAL::FE = FE_SimplexP<2>(3) +DEAL::using AMR +DEAL::L2 error = 0.003640 +DEAL::ratio = 274.7 +DEAL::L2 error = 0.0002153 +DEAL::ratio = 16.91 +DEAL::L2 error = 1.432e-05 +DEAL::ratio = 15.03 -- 2.39.5