From 8ac12522ba9256ba379e2732e60c65fc86358030 Mon Sep 17 00:00:00 2001 From: Marc Fehling Date: Mon, 25 Jan 2021 14:28:15 -0700 Subject: [PATCH] simplex: implement 'interface_constraints' for FE_Poly. --- source/simplex/fe_lib.cc | 67 +++++++++++++++++++ tests/simplex/hanging_nodes_01.cc | 2 +- ...ng_nodes_01.with_simplex_support=on.output | 2 + 3 files changed, 70 insertions(+), 1 deletion(-) diff --git a/source/simplex/fe_lib.cc b/source/simplex/fe_lib.cc index 61c14b5a43..585ab8affb 100644 --- a/source/simplex/fe_lib.cc +++ b/source/simplex/fe_lib.cc @@ -151,6 +151,72 @@ namespace Simplex return unit_face_points; } + /** + * Specify the constraints which the dofs on the two sides of a cell + * interface underlie if the line connects two cells of which one is refined + * once. + */ + template + FullMatrix + constraints_fe_poly(const unsigned int /*degree*/) + { + // no constraints in 1d + // constraints in 3d not implemented yet + return FullMatrix(); + } + + template <> + FullMatrix + constraints_fe_poly<2>(const unsigned int degree) + { + const unsigned int dim = 2; + + Assert(degree <= 2, ExcNotImplemented()); + + // the following implements the 2d case + // (the 3d case is not implemented yet) + // + // consult FE_Q_Base::Implementation::initialize_constraints() + // for more information + + 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); + } + + // Now construct relation between destination (child) and source (mother) + // dofs. + + const unsigned int n_dofs_constrained = constraint_points.size(); + unsigned int n_dofs_per_face = degree + 1; + FullMatrix interface_constraints(n_dofs_constrained, + n_dofs_per_face); + + const auto poly = + Simplex::BarycentricPolynomials::get_fe_p_basis(degree); + + for (unsigned int i = 0; i < n_dofs_constrained; ++i) + for (unsigned int j = 0; j < n_dofs_per_face; ++j) + { + interface_constraints(i, j) = + poly.compute_value(j, constraint_points[i]); + + // if the value is small up to round-off, then simply set it to zero + // to avoid unwanted fill-in of the constraint matrices (which would + // then increase the number of other DoFs a constrained DoF would + // couple to) + if (std::fabs(interface_constraints(i, j)) < 1e-13) + interface_constraints(i, j) = 0; + } + return interface_constraints; + } + /** * Helper function to set up the dpo vector of FE_DGP for a given @p dim and * @p degree. @@ -312,6 +378,7 @@ namespace Simplex if (conformity == FiniteElementData::Conformity::H1) this->unit_face_support_points = unit_face_support_points_fe_poly(degree); + this->interface_constraints = constraints_fe_poly(degree); } diff --git a/tests/simplex/hanging_nodes_01.cc b/tests/simplex/hanging_nodes_01.cc index 753ecdf45e..afbd51fbb9 100644 --- a/tests/simplex/hanging_nodes_01.cc +++ b/tests/simplex/hanging_nodes_01.cc @@ -53,7 +53,7 @@ test() // hanging node constraints AffineConstraints constraints; - // DoFTools::make_hanging_node_constraints(dofh, constraints); + DoFTools::make_hanging_node_constraints(dofh, constraints); constraints.print(deallog.get_file_stream()); deallog << "OK" << std::endl; diff --git a/tests/simplex/hanging_nodes_01.with_simplex_support=on.output b/tests/simplex/hanging_nodes_01.with_simplex_support=on.output index 1f03256372..8034a9e256 100644 --- a/tests/simplex/hanging_nodes_01.with_simplex_support=on.output +++ b/tests/simplex/hanging_nodes_01.with_simplex_support=on.output @@ -1,3 +1,5 @@ DEAL:2d::ndofs: 7 + 6 2: 0.500000 + 6 1: 0.500000 DEAL:2d::OK -- 2.39.5