From d36c8ba3b6a1712c6171aad9b8d8afee8492af4f Mon Sep 17 00:00:00 2001 From: Nuno Nobre Date: Tue, 23 Apr 2024 15:57:15 +0100 Subject: [PATCH] =?UTF-8?q?Provision=20for=20high-order=20N=C3=A9d=C3=A9le?= =?UTF-8?q?c=20elements?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/fe/fe_nedelec_one.C | 524 +++++++++++++++------------------------- 1 file changed, 194 insertions(+), 330 deletions(-) diff --git a/src/fe/fe_nedelec_one.C b/src/fe/fe_nedelec_one.C index 4b2c0a3bbad..6e26cc7bb58 100644 --- a/src/fe/fe_nedelec_one.C +++ b/src/fe/fe_nedelec_one.C @@ -39,6 +39,8 @@ LIBMESH_DEFAULT_VECTORIZED_FE(3,NEDELEC_ONE) // Anonymous namespace for local helper functions namespace { +// Forward-declare nedelec_one_n_dofs for immediate use +unsigned int nedelec_one_n_dofs(const ElemType, const Order); void nedelec_one_nodal_soln(const Elem * elem, const Order order, const std::vector & elem_soln, @@ -53,109 +55,43 @@ void nedelec_one_nodal_soln(const Elem * elem, nodal_soln.resize(n_nodes*dim); - FEType fe_type(order, NEDELEC_ONE); FEType p_refined_fe_type(totalorder, NEDELEC_ONE); - switch (totalorder) - { - case FIRST: - { - switch (elem_type) - { - case TRI6: - case TRI7: - { - libmesh_assert_equal_to (elem_soln.size(), 3); - - if (elem_type == TRI6) - libmesh_assert_equal_to (nodal_soln.size(), 6*2); - else - libmesh_assert_equal_to (nodal_soln.size(), 7*2); - break; - } - case QUAD8: - case QUAD9: - { - libmesh_assert_equal_to (elem_soln.size(), 4); - - if (elem_type == QUAD8) - libmesh_assert_equal_to (nodal_soln.size(), 8*2); - else - libmesh_assert_equal_to (nodal_soln.size(), 9*2); - break; - } - case TET10: - case TET14: - { - libmesh_assert_equal_to (elem_soln.size(), 6); - if (elem_type == TET10) - libmesh_assert_equal_to (nodal_soln.size(), 10*3); - else - libmesh_assert_equal_to (nodal_soln.size(), 14*3); - - break; - } - case HEX20: - case HEX27: - { - libmesh_assert_equal_to (elem_soln.size(), 12); - - if (elem_type == HEX20) - libmesh_assert_equal_to (nodal_soln.size(), 20*3); - else - libmesh_assert_equal_to (nodal_soln.size(), 27*3); - - break; - } - - default: - libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(elem_type) << " selected for NEDELEC_ONE FE family!"); - - } // switch(elem_type) + if (elem_type != TRI6 && elem_type != TRI7 && + elem_type != QUAD8 && elem_type != QUAD9 && + elem_type != TET10 && elem_type != TET14 && + elem_type != HEX20 && elem_type != HEX27) + libmesh_error_msg("ERROR: Invalid ElemType " << Utility::enum_to_string(elem_type) << " selected for NEDELEC_ONE FE family!"); - const unsigned int n_sf = - FEInterface::n_shape_functions(fe_type, elem); + libmesh_assert_equal_to (elem_soln.size(), nedelec_one_n_dofs(elem_type, totalorder)); - std::vector refspace_nodes; - FEVectorBase::get_refspace_nodes(elem_type,refspace_nodes); - libmesh_assert_equal_to (refspace_nodes.size(), n_nodes); + const unsigned int n_sf = FEInterface::n_shape_functions(p_refined_fe_type, elem); + std::vector refspace_nodes; + FEVectorBase::get_refspace_nodes(elem_type,refspace_nodes); + libmesh_assert_equal_to (refspace_nodes.size(), n_nodes); - // Need to create new fe object so the shape function has the FETransformation - // applied to it. - std::unique_ptr vis_fe = FEVectorBase::build(dim, p_refined_fe_type); + // Need to create new fe object so the shape function has the FETransformation + // applied to it. + std::unique_ptr vis_fe = FEVectorBase::build(dim, p_refined_fe_type); - const std::vector> & vis_phi = vis_fe->get_phi(); + const std::vector> & vis_phi = vis_fe->get_phi(); - vis_fe->reinit(elem,&refspace_nodes); + vis_fe->reinit(elem,&refspace_nodes); - for (unsigned int n = 0; n < n_nodes; n++) - { - libmesh_assert_equal_to (elem_soln.size(), n_sf); - - // Zero before summation - for (int d = 0; d < dim; d++) - { - nodal_soln[dim*n+d] = 0; - } - - // u = Sum (u_i phi_i) - for (unsigned int i=0; i unsigned int FE<1,NEDELEC_ONE>::n_dofs(const ElemType, const Order) template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); } template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs(const ElemType t, const Order o) { return nedelec_one_n_dofs(t, o); } - -// Do full-specialization for every dimension, instead -// of explicit instantiation at the end of this function. template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE } template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { NEDELEC_LOW_D_ERROR_MESSAGE } template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); } template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return nedelec_one_n_dofs_at_node(t, o, n); } - -// Nedelec first type elements have no dofs per element -// FIXME: Only for first order! template <> unsigned int FE<0,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE } template <> unsigned int FE<1,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { NEDELEC_LOW_D_ERROR_MESSAGE } -template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } -template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const ElemType, const Order) { return 0; } +template <> unsigned int FE<2,NEDELEC_ONE>::n_dofs_per_elem(const ElemType t, const Order o) { return nedelec_one_n_dofs_per_elem(t, o); } +template <> unsigned int FE<3,NEDELEC_ONE>::n_dofs_per_elem(const ElemType t, const Order o) { return nedelec_one_n_dofs_per_elem(t, o); } // Nedelec first type FEMs are always tangentially continuous template <> FEContinuity FE<0,NEDELEC_ONE>::get_continuity() const { NEDELEC_LOW_D_ERROR_MESSAGE }