Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
       3             : 
       4             : // This library is free software; you can redistribute it and/or
       5             : // modify it under the terms of the GNU Lesser General Public
       6             : // License as published by the Free Software Foundation; either
       7             : // version 2.1 of the License, or (at your option) any later version.
       8             : 
       9             : // This library is distributed in the hope that it will be useful,
      10             : // but WITHOUT ANY WARRANTY; without even the implied warranty of
      11             : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      12             : // Lesser General Public License for more details.
      13             : 
      14             : // You should have received a copy of the GNU Lesser General Public
      15             : // License along with this library; if not, write to the Free Software
      16             : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
      17             : 
      18             : 
      19             : 
      20             : // Local includes
      21             : #include "libmesh/libmesh_config.h"
      22             : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
      23             : 
      24             : #include "libmesh/elem.h"
      25             : #include "libmesh/enum_to_string.h"
      26             : #include "libmesh/fe.h"
      27             : #include "libmesh/fe_interface.h"
      28             : #include "libmesh/fe_macro.h"
      29             : 
      30             : // Anonymous namespace for local helper functions
      31             : namespace {
      32             : 
      33             : using namespace libMesh;
      34             : 
      35             : static const FEFamily _underlying_fe_family = BERNSTEIN;
      36             : 
      37      204160 : void rational_nodal_soln(const Elem * elem,
      38             :                          const Order order,
      39             :                          const std::vector<Number> & elem_soln,
      40             :                          std::vector<Number> & nodal_soln,
      41             :                          const bool add_p_level)
      42             : {
      43      204160 :   const unsigned int n_nodes = elem->n_nodes();
      44             : 
      45      204160 :   const ElemType elem_type = elem->type();
      46             : 
      47      204160 :   nodal_soln.resize(n_nodes);
      48             : 
      49      255200 :   const Order totalorder = order + add_p_level*elem->p_level();
      50             : 
      51             :   // FEType object to be passed to various FEInterface functions below.
      52       51040 :   FEType fe_type(order, _underlying_fe_family);
      53       51040 :   FEType p_refined_fe_type(totalorder, _underlying_fe_family);
      54             : 
      55      204160 :   switch (totalorder)
      56             :     {
      57             :       // Constant shape functions
      58           0 :     case CONSTANT:
      59             :       {
      60           0 :         libmesh_assert_equal_to (elem_soln.size(), 1);
      61             : 
      62           0 :         std::fill(nodal_soln.begin(), nodal_soln.end(), elem_soln[0]);
      63             : 
      64           0 :         return;
      65             :       }
      66             : 
      67             : 
      68             :       // For other bases do interpolation at the nodes
      69             :       // explicitly.
      70      204160 :     default:
      71             :       {
      72             :         const unsigned int n_sf =
      73      204160 :           FEInterface::n_shape_functions(fe_type, elem);
      74             : 
      75      102080 :         std::vector<Point> refspace_nodes;
      76      204160 :         FEBase::get_refspace_nodes(elem_type,refspace_nodes);
      77       51040 :         libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
      78       51040 :         libmesh_assert_equal_to (n_sf, n_nodes);
      79       51040 :         libmesh_assert_equal_to (elem_soln.size(), n_sf);
      80             : 
      81      255200 :         std::vector<Real> node_weights(n_nodes);
      82             : 
      83      102080 :         const unsigned char weight_index = elem->mapping_data();
      84             : 
      85     1129216 :         for (unsigned int n=0; n<n_nodes; n++)
      86     1156320 :           node_weights[n] =
      87      925056 :             elem->node_ref(n).get_extra_datum<Real>(weight_index);
      88             : 
      89             :         // Zero before summation
      90       51040 :         std::fill(nodal_soln.begin(), nodal_soln.end(), 0);
      91             : 
      92     1129216 :         for (unsigned int n=0; n<n_nodes; n++)
      93             :           {
      94      925056 :             std::vector<Real> weighted_shape(n_sf);
      95      231264 :             Real weighted_sum = 0;
      96             : 
      97     8338176 :             for (unsigned int i=0; i<n_sf; i++)
      98             :               {
      99     7413120 :                 weighted_shape[i] = node_weights[i] *
     100     9266400 :                   FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
     101     7413120 :                 weighted_sum += weighted_shape[i];
     102             :               }
     103             : 
     104             :             // u_i = Sum (alpha_i w_i phi_i) / Sum (w_j phi_j)
     105     8338176 :             for (unsigned int i=0; i<n_sf; i++)
     106    12972960 :               nodal_soln[n] += elem_soln[i] * weighted_shape[i];
     107     1156320 :             nodal_soln[n] /= weighted_sum;
     108             :           }
     109             : 
     110       51040 :         return;
     111             :       }
     112             :     }
     113             : } //  rational_nodal_soln()
     114             : 
     115             : } // anon namespace
     116             : 
     117             : 
     118             : namespace libMesh {
     119             : 
     120             : 
     121             : // Instantiate (side_) nodal_soln() function for every dimension
     122      204160 : LIBMESH_FE_NODAL_SOLN(RATIONAL_BERNSTEIN, rational_nodal_soln)
     123           0 : LIBMESH_FE_SIDE_NODAL_SOLN(RATIONAL_BERNSTEIN)
     124             : 
     125             : 
     126             : // Full specialization of n_dofs() function for every dimension
     127           0 : template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return FE<0,_underlying_fe_family>::n_dofs(t, o); }
     128           0 : template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return FE<1,_underlying_fe_family>::n_dofs(t, o); }
     129           0 : template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return FE<2,_underlying_fe_family>::n_dofs(t, o); }
     130           0 : template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs(const ElemType t, const Order o) { return FE<3,_underlying_fe_family>::n_dofs(t, o); }
     131             : 
     132      602469 : template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<0,_underlying_fe_family>::n_dofs(e, o); }
     133      100895 : template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<1,_underlying_fe_family>::n_dofs(e, o); }
     134     6075827 : template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<2,_underlying_fe_family>::n_dofs(e, o); }
     135     1097953 : template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs(const Elem * e, const Order o) { return FE<3,_underlying_fe_family>::n_dofs(e, o); }
     136             : 
     137             : // Full specialization of n_dofs_at_node() function for every dimension.
     138     1635220 : template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<0,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
     139        3798 : template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<1,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
     140    12100167 : template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<2,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
     141      559107 : template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return FE<3,_underlying_fe_family>::n_dofs_at_node(t, o, n); }
     142             : 
     143      126200 : template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<0,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
     144        1647 : template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<1,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
     145      597758 : template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<2,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
     146      290514 : template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_at_node(const Elem & e, const Order o, const unsigned int n) { return FE<3,_underlying_fe_family>::n_dofs_at_node(e, o, n); }
     147             : 
     148             : // Full specialization of n_dofs_per_elem() function for every dimension.
     149           0 : template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<0,_underlying_fe_family>::n_dofs_per_elem(t, o); }
     150           0 : template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<1,_underlying_fe_family>::n_dofs_per_elem(t, o); }
     151           0 : template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<2,_underlying_fe_family>::n_dofs_per_elem(t, o); }
     152           0 : template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const ElemType t, const Order o) { return FE<3,_underlying_fe_family>::n_dofs_per_elem(t, o); }
     153             : 
     154     1698320 : template <> unsigned int FE<0,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<0,_underlying_fe_family>::n_dofs_per_elem(e, o); }
     155        1962 : template <> unsigned int FE<1,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<1,_underlying_fe_family>::n_dofs_per_elem(e, o); }
     156     1386721 : template <> unsigned int FE<2,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<2,_underlying_fe_family>::n_dofs_per_elem(e, o); }
     157       53439 : template <> unsigned int FE<3,RATIONAL_BERNSTEIN>::n_dofs_per_elem(const Elem & e, const Order o) { return FE<3,_underlying_fe_family>::n_dofs_per_elem(e, o); }
     158             : 
     159             : // Our current rational FEMs are C^0 continuous
     160        1618 : template <> FEContinuity FE<0,RATIONAL_BERNSTEIN>::get_continuity() const { return C_ZERO; }
     161        6576 : template <> FEContinuity FE<1,RATIONAL_BERNSTEIN>::get_continuity() const { return C_ZERO; }
     162       74357 : template <> FEContinuity FE<2,RATIONAL_BERNSTEIN>::get_continuity() const { return C_ZERO; }
     163       18961 : template <> FEContinuity FE<3,RATIONAL_BERNSTEIN>::get_continuity() const { return C_ZERO; }
     164             : 
     165             : // Rational FEMs are in general not hierarchic
     166           0 : template <> bool FE<0,RATIONAL_BERNSTEIN>::is_hierarchic() const { return false; }
     167          24 : template <> bool FE<1,RATIONAL_BERNSTEIN>::is_hierarchic() const { return false; }
     168          46 : template <> bool FE<2,RATIONAL_BERNSTEIN>::is_hierarchic() const { return false; }
     169          90 : template <> bool FE<3,RATIONAL_BERNSTEIN>::is_hierarchic() const { return false; }
     170             : 
     171             : #ifdef LIBMESH_ENABLE_AMR
     172             : // compute_constraints() specializations are only needed for 2 and 3D
     173             : template <>
     174       52726 : void FE<2,RATIONAL_BERNSTEIN>::compute_constraints (DofConstraints & constraints,
     175             :                                                       DofMap & dof_map,
     176             :                                                       const unsigned int variable_number,
     177             :                                                       const Elem * elem)
     178       52726 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
     179             : 
     180             : template <>
     181        3016 : void FE<3,RATIONAL_BERNSTEIN>::compute_constraints (DofConstraints & constraints,
     182             :                                                       DofMap & dof_map,
     183             :                                                       const unsigned int variable_number,
     184             :                                                       const Elem * elem)
     185        3016 : { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
     186             : #endif // #ifdef LIBMESH_ENABLE_AMR
     187             : 
     188             : // Bernstein shapes need reinit only for approximation orders >= 3,
     189             : // but for *RATIONAL* shapes we could need reinit anywhere because our
     190             : // nodal weights can change from element to element.
     191      195696 : template <> bool FE<0,RATIONAL_BERNSTEIN>::shapes_need_reinit() const { return true; }
     192         180 : template <> bool FE<1,RATIONAL_BERNSTEIN>::shapes_need_reinit() const { return true; }
     193      216521 : template <> bool FE<2,RATIONAL_BERNSTEIN>::shapes_need_reinit() const { return true; }
     194        6184 : template <> bool FE<3,RATIONAL_BERNSTEIN>::shapes_need_reinit() const { return true; }
     195             : 
     196             : } // namespace libMesh
     197             : 
     198             : #endif //LIBMESH_ENABLE_HIGHER_ORDER_SHAPES