1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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/elem.h"
22 #include "libmesh/fe.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/enum_to_string.h"
25 
26 namespace libMesh
27 {
28 
29 // ------------------------------------------------------------
30 // Clough-specific implementations
31 
32 // Anonymous namespace for local helper functions
33 namespace {
34 
clough_nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)35 void clough_nodal_soln(const Elem * elem,
36                        const Order order,
37                        const std::vector<Number> & elem_soln,
38                        std::vector<Number> &       nodal_soln)
39 {
40   const unsigned int n_nodes = elem->n_nodes();
41 
42   const ElemType elem_type = elem->type();
43 
44   nodal_soln.resize(n_nodes);
45 
46   const Order totalorder = static_cast<Order>(order + elem->p_level());
47 
48   // FEType object to be passed to various FEInterface functions below.
49   FEType fe_type(order, CLOUGH);
50 
51   switch (totalorder)
52     {
53       // Piecewise cubic shape functions with linear flux edges
54     case SECOND:
55       // Piecewise cubic shape functions
56     case THIRD:
57       {
58         const unsigned int n_sf =
59           FEInterface::n_shape_functions(fe_type, elem);
60 
61         std::vector<Point> refspace_nodes;
62         FEBase::get_refspace_nodes(elem_type,refspace_nodes);
63         libmesh_assert_equal_to (refspace_nodes.size(), n_nodes);
64 
65         for (unsigned int n=0; n<n_nodes; n++)
66           {
67             libmesh_assert_equal_to (elem_soln.size(), n_sf);
68 
69             // Zero before summation
70             nodal_soln[n] = 0;
71 
72             // u_i = Sum (alpha_i phi_i)
73             for (unsigned int i=0; i<n_sf; i++)
74               nodal_soln[n] += elem_soln[i] *
75                 FEInterface::shape(fe_type, elem, i, refspace_nodes[n]);
76           }
77 
78         return;
79       }
80 
81     default:
82       libmesh_error_msg("ERROR: Invalid total order " << totalorder);
83     }
84 } // clough_nodal_soln()
85 
86 
87 
88 
clough_n_dofs(const ElemType t,const Order o)89 unsigned int clough_n_dofs(const ElemType t, const Order o)
90 {
91   if (t == INVALID_ELEM)
92     return 0;
93 
94   switch (o)
95     {
96       // Piecewise cubic shape functions with linear flux edges
97     case SECOND:
98       {
99         switch (t)
100           {
101           case TRI6:
102             return 9;
103 
104           default:
105             libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
106           }
107       }
108       // Piecewise cubic Clough-Tocher element
109     case THIRD:
110       {
111         switch (t)
112           {
113           case TRI6:
114             return 12;
115 
116           default:
117             libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
118           }
119       }
120 
121     default:
122       libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
123     }
124 } // clough_n_dofs()
125 
126 
127 
128 
129 
clough_n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)130 unsigned int clough_n_dofs_at_node(const ElemType t,
131                                    const Order o,
132                                    const unsigned int n)
133 {
134   switch (o)
135     {
136       // Piecewise cubic shape functions with linear flux edges
137     case SECOND:
138       {
139         switch (t)
140           {
141             // The 2D Clough-Tocher defined on a 6-noded triangle
142           case TRI6:
143             {
144               switch (n)
145                 {
146                 case 0:
147                 case 1:
148                 case 2:
149                   return 3;
150 
151                 case 3:
152                 case 4:
153                 case 5:
154                   return 0;
155 
156                 default:
157                   libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
158                 }
159             }
160 
161           default:
162             libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
163           }
164       }
165       // The third-order clough shape functions
166     case THIRD:
167       {
168         switch (t)
169           {
170             // The 2D Clough-Tocher defined on a 6-noded triangle
171           case TRI6:
172             {
173               switch (n)
174                 {
175                 case 0:
176                 case 1:
177                 case 2:
178                   return 3;
179 
180                 case 3:
181                 case 4:
182                 case 5:
183                   return 1;
184 
185                 default:
186                   libmesh_error_msg("ERROR: Invalid node ID " << n << " selected for TRI6!");
187                 }
188             }
189 
190           default:
191             libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
192           }
193       }
194     default:
195       libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
196     }
197 } // clough_n_dofs_at_node()
198 
199 
200 
201 
clough_n_dofs_per_elem(const ElemType t,const Order o)202 unsigned int clough_n_dofs_per_elem(const ElemType t, const Order o)
203 {
204   switch (o)
205     {
206       // Piecewise cubic shape functions with linear flux edges
207     case SECOND:
208       // The third-order Clough-Tocher shape functions
209     case THIRD:
210       {
211         switch (t)
212           {
213             // The 2D clough defined on a 6-noded triangle
214           case TRI6:
215             return 0;
216 
217           default:
218             libmesh_error_msg("ERROR: Bad ElemType = " << Utility::enum_to_string(t) << " for " << Utility::enum_to_string(o) << " order approximation!");
219           }
220       }
221 
222     default:
223       libmesh_error_msg("ERROR: Invalid Order " << Utility::enum_to_string(o) << " selected for CLOUGH FE family!");
224     }
225 } // clough_n_dofs_per_elem()
226 
227 } // anonymous
228 
229 
230 
231 
232   // Do full-specialization of nodal_soln() function for every
233   // dimension, instead of explicit instantiation at the end of this
234   // file.
235   // This could be macro-ified so that it fits on one line...
236 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)237 void FE<0,CLOUGH>::nodal_soln(const Elem * elem,
238                               const Order order,
239                               const std::vector<Number> & elem_soln,
240                               std::vector<Number> & nodal_soln)
241 { clough_nodal_soln(elem, order, elem_soln, nodal_soln); }
242 
243 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)244 void FE<1,CLOUGH>::nodal_soln(const Elem * elem,
245                               const Order order,
246                               const std::vector<Number> & elem_soln,
247                               std::vector<Number> & nodal_soln)
248 { clough_nodal_soln(elem, order, elem_soln, nodal_soln); }
249 
250 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)251 void FE<2,CLOUGH>::nodal_soln(const Elem * elem,
252                               const Order order,
253                               const std::vector<Number> & elem_soln,
254                               std::vector<Number> & nodal_soln)
255 { clough_nodal_soln(elem, order, elem_soln, nodal_soln); }
256 
257 template <>
nodal_soln(const Elem * elem,const Order order,const std::vector<Number> & elem_soln,std::vector<Number> & nodal_soln)258 void FE<3,CLOUGH>::nodal_soln(const Elem * elem,
259                               const Order order,
260                               const std::vector<Number> & elem_soln,
261                               std::vector<Number> & nodal_soln)
262 { clough_nodal_soln(elem, order, elem_soln, nodal_soln); }
263 
264 
265 // Full specialization of n_dofs() function for every dimension
n_dofs(const ElemType t,const Order o)266 template <> unsigned int FE<0,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)267 template <> unsigned int FE<1,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)268 template <> unsigned int FE<2,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
n_dofs(const ElemType t,const Order o)269 template <> unsigned int FE<3,CLOUGH>::n_dofs(const ElemType t, const Order o) { return clough_n_dofs(t, o); }
270 
271 
272 // Full specialization of n_dofs_at_node() function for every dimension.
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)273 template <> unsigned int FE<0,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)274 template <> unsigned int FE<1,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)275 template <> unsigned int FE<2,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
n_dofs_at_node(const ElemType t,const Order o,const unsigned int n)276 template <> unsigned int FE<3,CLOUGH>::n_dofs_at_node(const ElemType t, const Order o, const unsigned int n) { return clough_n_dofs_at_node(t, o, n); }
277 
278 // Full specialization of n_dofs_per_elem() function for every dimension.
n_dofs_per_elem(const ElemType t,const Order o)279 template <> unsigned int FE<0,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
n_dofs_per_elem(const ElemType t,const Order o)280 template <> unsigned int FE<1,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
n_dofs_per_elem(const ElemType t,const Order o)281 template <> unsigned int FE<2,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
n_dofs_per_elem(const ElemType t,const Order o)282 template <> unsigned int FE<3,CLOUGH>::n_dofs_per_elem(const ElemType t, const Order o) { return clough_n_dofs_per_elem(t, o); }
283 
284 // Clough FEMs are C^1 continuous
get_continuity()285 template <> FEContinuity FE<0,CLOUGH>::get_continuity() const { return C_ONE; }
get_continuity()286 template <> FEContinuity FE<1,CLOUGH>::get_continuity() const { return C_ONE; }
get_continuity()287 template <> FEContinuity FE<2,CLOUGH>::get_continuity() const { return C_ONE; }
get_continuity()288 template <> FEContinuity FE<3,CLOUGH>::get_continuity() const { return C_ONE; }
289 
290 // Clough FEMs are not (currently) hierarchic
is_hierarchic()291 template <> bool FE<0,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
is_hierarchic()292 template <> bool FE<1,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
is_hierarchic()293 template <> bool FE<2,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
is_hierarchic()294 template <> bool FE<3,CLOUGH>::is_hierarchic() const { return false; } // FIXME - this will be changed
295 
296 #ifdef LIBMESH_ENABLE_AMR
297 // compute_constraints() specializations are only needed for 2 and 3D
298 template <>
compute_constraints(DofConstraints & constraints,DofMap & dof_map,const unsigned int variable_number,const Elem * elem)299 void FE<2,CLOUGH>::compute_constraints (DofConstraints & constraints,
300                                         DofMap & dof_map,
301                                         const unsigned int variable_number,
302                                         const Elem * elem)
303 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
304 
305 template <>
compute_constraints(DofConstraints & constraints,DofMap & dof_map,const unsigned int variable_number,const Elem * elem)306 void FE<3,CLOUGH>::compute_constraints (DofConstraints & constraints,
307                                         DofMap & dof_map,
308                                         const unsigned int variable_number,
309                                         const Elem * elem)
310 { compute_proj_constraints(constraints, dof_map, variable_number, elem); }
311 #endif // #ifdef LIBMESH_ENABLE_AMR
312 
313 // Clough FEM shapes need reinit
shapes_need_reinit()314 template <> bool FE<0,CLOUGH>::shapes_need_reinit() const { return true; }
shapes_need_reinit()315 template <> bool FE<1,CLOUGH>::shapes_need_reinit() const { return true; }
shapes_need_reinit()316 template <> bool FE<2,CLOUGH>::shapes_need_reinit() const { return true; }
shapes_need_reinit()317 template <> bool FE<3,CLOUGH>::shapes_need_reinit() const { return true; }
318 
319 } // namespace libMesh
320