diff --git a/include/geom/cell_tet14.h b/include/geom/cell_tet14.h index ac655270a40..ba615e103df 100644 --- a/include/geom/cell_tet14.h +++ b/include/geom/cell_tet14.h @@ -91,7 +91,7 @@ class Tet14 final : public Tet virtual ElemType type () const override { return TET14; } /** - * \returns 10. + * \returns 14. */ virtual unsigned int n_nodes() const override { return num_nodes; } diff --git a/src/fe/fe_hierarchic.C b/src/fe/fe_hierarchic.C index 922d1d49ea0..ca5d2d71f10 100644 --- a/src/fe/fe_hierarchic.C +++ b/src/fe/fe_hierarchic.C @@ -117,8 +117,6 @@ unsigned int hierarchic_n_dofs(const ElemType t, const Order o) case QUAD9: return ((o+1)*(o+1)); case HEX8: - libmesh_assert_less (o, 2); - libmesh_fallthrough(); case HEX20: libmesh_assert_less (o, 2); libmesh_fallthrough(); diff --git a/src/fe/fe_nedelec_one.C b/src/fe/fe_nedelec_one.C index 4b2c0a3bbad..ff6b42c0a18 100644 --- a/src/fe/fe_nedelec_one.C +++ b/src/fe/fe_nedelec_one.C @@ -39,6 +39,10 @@ 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 +57,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 } diff --git a/src/fe/fe_nedelec_one_shape_2D.C b/src/fe/fe_nedelec_one_shape_2D.C index 4a9249275ce..fdb8ac999ad 100644 --- a/src/fe/fe_nedelec_one_shape_2D.C +++ b/src/fe/fe_nedelec_one_shape_2D.C @@ -37,6 +37,10 @@ RealGradient FE<2,NEDELEC_ONE>::shape(const Elem * elem, libmesh_assert(elem); const Order total_order = static_cast(order + add_p_level * elem->p_level()); + const short sign = elem->point(i) > elem->point((i+1) % elem->n_vertices()) ? 1 : -1; + + const Real xi = p(0); + const Real eta = p(1); switch (total_order) { @@ -50,9 +54,6 @@ RealGradient FE<2,NEDELEC_ONE>::shape(const Elem * elem, { libmesh_assert_less (i, 4); - const Real xi = p(0); - const Real eta = p(1); - // Even with a loose inverse_map tolerance we ought to // be nearly on the element interior in master // coordinates @@ -62,34 +63,13 @@ RealGradient FE<2,NEDELEC_ONE>::shape(const Elem * elem, switch(i) { case 0: - { - if (elem->point(0) > elem->point(1)) - return RealGradient( -0.25*(1.0-eta), 0.0 ); - else - return RealGradient( 0.25*(1.0-eta), 0.0 ); - } + return sign * RealGradient( -0.25*(1.0-eta), 0.0 ); case 1: - { - if (elem->point(1) > elem->point(2)) - return RealGradient( 0.0, -0.25*(1.0+xi) ); - else - return RealGradient( 0.0, 0.25*(1.0+xi) ); - } - + return sign * RealGradient( 0.0, -0.25*(1.0+xi) ); case 2: - { - if (elem->point(2) > elem->point(3)) - return RealGradient( 0.25*(1.0+eta), 0.0 ); - else - return RealGradient( -0.25*(1.0+eta), 0.0 ); - } + return sign * RealGradient( 0.25*(1.0+eta), 0.0 ); case 3: - { - if (elem->point(3) > elem->point(0)) - return RealGradient( 0.0, -0.25*(xi-1.0) ); - else - return RealGradient( 0.0, 0.25*(xi-1.0) ); - } + return sign * RealGradient( 0.0, -0.25*(xi-1.0) ); default: libmesh_error_msg("Invalid i = " << i); @@ -101,35 +81,16 @@ RealGradient FE<2,NEDELEC_ONE>::shape(const Elem * elem, case TRI6: case TRI7: { - const Real xi = p(0); - const Real eta = p(1); - libmesh_assert_less (i, 3); switch(i) { case 0: - { - if (elem->point(0) > elem->point(1)) - return RealGradient( -1.0+eta, -xi ); - else - return RealGradient( 1.0-eta, xi ); - } + return sign * RealGradient( -1.0+eta, -xi ); case 1: - { - if (elem->point(1) > elem->point(2)) - return RealGradient( eta, -xi ); - else - return RealGradient( -eta, xi ); - } - + return sign * RealGradient( eta, -xi ); case 2: - { - if (elem->point(2) > elem->point(0)) - return RealGradient( eta, -xi+1.0 ); - else - return RealGradient( -eta, xi-1.0 ); - } + return sign * RealGradient( eta, -xi+1.0 ); default: libmesh_error_msg("Invalid i = " << i); @@ -188,6 +149,7 @@ RealGradient FE<2,NEDELEC_ONE>::shape_deriv(const Elem * elem, libmesh_assert_less (j, 2); const Order total_order = static_cast(order + add_p_level * elem->p_level()); + const short sign = elem->point(i) > elem->point((i+1) % elem->n_vertices()) ? 1 : -1; switch (total_order) { @@ -212,19 +174,9 @@ RealGradient FE<2,NEDELEC_ONE>::shape_deriv(const Elem * elem, case 2: return RealGradient(); case 1: - { - if (elem->point(1) > elem->point(2)) - return RealGradient( 0.0, -0.25 ); - else - return RealGradient( 0.0, 0.25 ); - } case 3: - { - if (elem->point(3) > elem->point(0)) - return RealGradient( 0.0, -0.25 ); - else - return RealGradient( 0.0, 0.25 ); - } + return sign * RealGradient( 0.0, -0.25 ); + default: libmesh_error_msg("Invalid i = " << i); } @@ -239,19 +191,9 @@ RealGradient FE<2,NEDELEC_ONE>::shape_deriv(const Elem * elem, case 3: return RealGradient(); case 0: - { - if (elem->point(0) > elem->point(1)) - return RealGradient( 0.25 ); - else - return RealGradient( -0.25 ); - } case 2: - { - if (elem->point(2) > elem->point(3)) - return RealGradient( 0.25 ); - else - return RealGradient( -0.25 ); - } + return sign * RealGradient( 0.25 ); + default: libmesh_error_msg("Invalid i = " << i); } @@ -260,8 +202,6 @@ RealGradient FE<2,NEDELEC_ONE>::shape_deriv(const Elem * elem, default: libmesh_error_msg("Invalid j = " << j); } - - return RealGradient(); } case TRI6: @@ -269,45 +209,16 @@ RealGradient FE<2,NEDELEC_ONE>::shape_deriv(const Elem * elem, { libmesh_assert_less (i, 3); - // Account for edge flipping - Real f = 1.0; - - switch(i) - { - case 0: - { - if (elem->point(0) > elem->point(1)) - f = -1.0; - break; - } - case 1: - { - if (elem->point(1) > elem->point(2)) - f = -1.0; - break; - } - case 2: - { - if (elem->point(2) > elem->point(0)) - f = -1.0; - break; - } - default: - libmesh_error_msg("Invalid i = " << i); - } - switch (j) { // d()/dxi case 0: - { - return RealGradient( 0.0, f*1.0); - } + return sign * RealGradient( 0.0, -1.0 ); + // d()/deta case 1: - { - return RealGradient( f*(-1.0) ); - } + return sign * RealGradient( 1.0 ); + default: libmesh_error_msg("Invalid j = " << j); }