fe_l2_lagrange.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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/dof_map.h"
22 #include "libmesh/fe.h"
23 #include "libmesh/fe_interface.h"
24 #include "libmesh/elem.h"
25 #include "libmesh/threads.h"
26 
27 namespace libMesh
28 {
29 
30  // ------------------------------------------------------------
31  // Lagrange-specific implementations
32 
33 
34  // Anonymous namespace for local helper functions
35  namespace {
36  void l2_lagrange_nodal_soln(const Elem* elem,
37  const Order order,
38  const std::vector<Number>& elem_soln,
39  std::vector<Number>& nodal_soln)
40  {
41  const unsigned int n_nodes = elem->n_nodes();
42  const ElemType type = elem->type();
43 
44  const Order totalorder = static_cast<Order>(order+elem->p_level());
45 
46  nodal_soln.resize(n_nodes);
47 
48 
49 
50  switch (totalorder)
51  {
52  // linear Lagrange shape functions
53  case FIRST:
54  {
55  switch (type)
56  {
57  case EDGE3:
58  {
59  libmesh_assert_equal_to (elem_soln.size(), 2);
60  libmesh_assert_equal_to (nodal_soln.size(), 3);
61 
62  nodal_soln[0] = elem_soln[0];
63  nodal_soln[1] = elem_soln[1];
64  nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]);
65 
66  return;
67  }
68 
69  case EDGE4:
70  {
71  libmesh_assert_equal_to (elem_soln.size(), 2);
72  libmesh_assert_equal_to (nodal_soln.size(), 4);
73 
74  nodal_soln[0] = elem_soln[0];
75  nodal_soln[1] = elem_soln[1];
76  nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.;
77  nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.;
78 
79  return;
80  }
81 
82 
83  case TRI6:
84  {
85  libmesh_assert_equal_to (elem_soln.size(), 3);
86  libmesh_assert_equal_to (nodal_soln.size(), 6);
87 
88  nodal_soln[0] = elem_soln[0];
89  nodal_soln[1] = elem_soln[1];
90  nodal_soln[2] = elem_soln[2];
91  nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]);
92  nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]);
93  nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]);
94 
95  return;
96  }
97 
98 
99  case QUAD8:
100  case QUAD9:
101  {
102  libmesh_assert_equal_to (elem_soln.size(), 4);
103 
104  if (type == QUAD8)
105  libmesh_assert_equal_to (nodal_soln.size(), 8);
106  else
107  libmesh_assert_equal_to (nodal_soln.size(), 9);
108 
109 
110  nodal_soln[0] = elem_soln[0];
111  nodal_soln[1] = elem_soln[1];
112  nodal_soln[2] = elem_soln[2];
113  nodal_soln[3] = elem_soln[3];
114  nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
115  nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
116  nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]);
117  nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
118 
119  if (type == QUAD9)
120  nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
121 
122  return;
123  }
124 
125 
126  case TET10:
127  {
128  libmesh_assert_equal_to (elem_soln.size(), 4);
129  libmesh_assert_equal_to (nodal_soln.size(), 10);
130 
131  nodal_soln[0] = elem_soln[0];
132  nodal_soln[1] = elem_soln[1];
133  nodal_soln[2] = elem_soln[2];
134  nodal_soln[3] = elem_soln[3];
135  nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
136  nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
137  nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]);
138  nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
139  nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]);
140  nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]);
141 
142  return;
143  }
144 
145 
146  case HEX20:
147  case HEX27:
148  {
149  libmesh_assert_equal_to (elem_soln.size(), 8);
150 
151  if (type == HEX20)
152  libmesh_assert_equal_to (nodal_soln.size(), 20);
153  else
154  libmesh_assert_equal_to (nodal_soln.size(), 27);
155 
156  nodal_soln[0] = elem_soln[0];
157  nodal_soln[1] = elem_soln[1];
158  nodal_soln[2] = elem_soln[2];
159  nodal_soln[3] = elem_soln[3];
160  nodal_soln[4] = elem_soln[4];
161  nodal_soln[5] = elem_soln[5];
162  nodal_soln[6] = elem_soln[6];
163  nodal_soln[7] = elem_soln[7];
164  nodal_soln[8] = .5*(elem_soln[0] + elem_soln[1]);
165  nodal_soln[9] = .5*(elem_soln[1] + elem_soln[2]);
166  nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]);
167  nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]);
168  nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]);
169  nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]);
170  nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]);
171  nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]);
172  nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]);
173  nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]);
174  nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]);
175  nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]);
176 
177  if (type == HEX27)
178  {
179  nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
180  nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]);
181  nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]);
182  nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]);
183  nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]);
184  nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
185 
186  nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] +
187  elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
188  }
189 
190  return;
191  }
192 
193 
194  case PRISM15:
195  case PRISM18:
196  {
197  libmesh_assert_equal_to (elem_soln.size(), 6);
198 
199  if (type == PRISM15)
200  libmesh_assert_equal_to (nodal_soln.size(), 15);
201  else
202  libmesh_assert_equal_to (nodal_soln.size(), 18);
203 
204  nodal_soln[0] = elem_soln[0];
205  nodal_soln[1] = elem_soln[1];
206  nodal_soln[2] = elem_soln[2];
207  nodal_soln[3] = elem_soln[3];
208  nodal_soln[4] = elem_soln[4];
209  nodal_soln[5] = elem_soln[5];
210  nodal_soln[6] = .5*(elem_soln[0] + elem_soln[1]);
211  nodal_soln[7] = .5*(elem_soln[1] + elem_soln[2]);
212  nodal_soln[8] = .5*(elem_soln[0] + elem_soln[2]);
213  nodal_soln[9] = .5*(elem_soln[0] + elem_soln[3]);
214  nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
215  nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]);
216  nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
217  nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]);
218  nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]);
219 
220  if (type == PRISM18)
221  {
222  nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]);
223  nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]);
224  nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]);
225  }
226 
227  return;
228  }
229 
230 
231 
232  default:
233  {
234  // By default the element solution _is_ nodal,
235  // so just copy it.
236  nodal_soln = elem_soln;
237 
238  return;
239  }
240  }
241  }
242 
243  case SECOND:
244  {
245  switch (type)
246  {
247  case EDGE4:
248  {
249  libmesh_assert_equal_to (elem_soln.size(), 3);
250  libmesh_assert_equal_to (nodal_soln.size(), 4);
251 
252  // Project quadratic solution onto cubic element nodes
253  nodal_soln[0] = elem_soln[0];
254  nodal_soln[1] = elem_soln[1];
255  nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] +
256  8.*elem_soln[2])/9.;
257  nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] +
258  8.*elem_soln[2])/9.;
259  return;
260  }
261 
262  default:
263  {
264  // By default the element solution _is_ nodal,
265  // so just copy it.
266  nodal_soln = elem_soln;
267 
268  return;
269  }
270  }
271  }
272 
273 
274 
275 
276  default:
277  {
278  // By default the element solution _is_ nodal,
279  // so just copy it.
280  nodal_soln = elem_soln;
281 
282  return;
283  }
284  }
285  }
286 
287 
288  // TODO: We should make this work, for example, for SECOND on a TRI3
289  // (this is valid with L2_LAGRANGE, but not with LAGRANGE)
290  unsigned int l2_lagrange_n_dofs(const ElemType t, const Order o)
291  {
292  switch (o)
293  {
294 
295  // linear Lagrange shape functions
296  case FIRST:
297  {
298  switch (t)
299  {
300  case NODEELEM:
301  return 1;
302 
303  case EDGE2:
304  case EDGE3:
305  case EDGE4:
306  return 2;
307 
308  case TRI3:
309  case TRI6:
310  return 3;
311 
312  case QUAD4:
313  case QUAD8:
314  case QUAD9:
315  return 4;
316 
317  case TET4:
318  case TET10:
319  return 4;
320 
321  case HEX8:
322  case HEX20:
323  case HEX27:
324  return 8;
325 
326  case PRISM6:
327  case PRISM15:
328  case PRISM18:
329  return 6;
330 
331  case PYRAMID5:
332  case PYRAMID14:
333  return 5;
334 
335  default:
336  {
337 #ifdef DEBUG
338  libMesh::err << "ERROR: Bad ElemType = " << t
339  << " for FIRST order approximation!"
340  << std::endl;
341 #endif
342  libmesh_error();
343  }
344  }
345  }
346 
347 
348  // quadratic Lagrange shape functions
349  case SECOND:
350  {
351  switch (t)
352  {
353  case NODEELEM:
354  return 1;
355 
356  case EDGE3:
357  return 3;
358 
359  case TRI6:
360  return 6;
361 
362  case QUAD8:
363  return 8;
364 
365  case QUAD9:
366  return 9;
367 
368  case TET10:
369  return 10;
370 
371  case HEX20:
372  return 20;
373 
374  case HEX27:
375  return 27;
376 
377  case PRISM15:
378  return 15;
379 
380  case PRISM18:
381  return 18;
382 
383  default:
384  {
385 #ifdef DEBUG
386  libMesh::err << "ERROR: Bad ElemType = " << t
387  << " for SECOND order approximation!"
388  << std::endl;
389 #endif
390  libmesh_error();
391  }
392  }
393  }
394 
395  case THIRD:
396  {
397  switch (t)
398  {
399  case NODEELEM:
400  return 1;
401 
402  case EDGE4:
403  return 4;
404 
405  default:
406  {
407 #ifdef DEBUG
408  libMesh::err << "ERROR: Bad ElemType = " << t
409  << " for THIRD order approximation!"
410  << std::endl;
411 #endif
412  libmesh_error();
413  }
414  }
415  }
416 
417  default:
418  libmesh_error();
419  }
420 
421  libmesh_error();
422  return 0;
423  }
424 
425  } // anonymous namespace
426 
427 
428  // Do full-specialization for every dimension, instead
429  // of explicit instantiation at the end of this file.
430  // This could be macro-ified so that it fits on one line...
431  template <>
433  const Order order,
434  const std::vector<Number>& elem_soln,
435  std::vector<Number>& nodal_soln)
436  { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
437 
438  template <>
440  const Order order,
441  const std::vector<Number>& elem_soln,
442  std::vector<Number>& nodal_soln)
443  { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
444 
445  template <>
447  const Order order,
448  const std::vector<Number>& elem_soln,
449  std::vector<Number>& nodal_soln)
450  { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
451 
452  template <>
454  const Order order,
455  const std::vector<Number>& elem_soln,
456  std::vector<Number>& nodal_soln)
457  { l2_lagrange_nodal_soln(elem, order, elem_soln, nodal_soln); }
458 
459 
460  // Do full-specialization for every dimension, instead
461  // of explicit instantiation at the end of this function.
462  // This could be macro-ified.
463  template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
464  template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
465  template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
466  template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
467 
468 
469  // L2 Lagrange elements have all dofs on elements (hence none on nodes)
470  template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
471  template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
472  template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
473  template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 0; }
474 
475 
476  // L2 Lagrange elements have all dofs on elements
477  template <> unsigned int FE<0,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
478  template <> unsigned int FE<1,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
479  template <> unsigned int FE<2,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
480  template <> unsigned int FE<3,L2_LAGRANGE>::n_dofs_per_elem(const ElemType t, const Order o) { return l2_lagrange_n_dofs(t, o); }
481 
482  // L2 Lagrange FEMs are DISCONTINUOUS
487 
488  // Lagrange FEMs are not hierarchic
489  template <> bool FE<0,L2_LAGRANGE>::is_hierarchic() const { return false; }
490  template <> bool FE<1,L2_LAGRANGE>::is_hierarchic() const { return false; }
491  template <> bool FE<2,L2_LAGRANGE>::is_hierarchic() const { return false; }
492  template <> bool FE<3,L2_LAGRANGE>::is_hierarchic() const { return false; }
493 
494  // Lagrange FEM shapes do not need reinit (is this always true?)
495  template <> bool FE<0,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
496  template <> bool FE<1,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
497  template <> bool FE<2,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
498  template <> bool FE<3,L2_LAGRANGE>::shapes_need_reinit() const { return false; }
499 
500  // We don't need any constraints for this DISCONTINUOUS basis, hence these
501  // are NOOPs
502 #ifdef LIBMESH_ENABLE_AMR
503  template <>
505  DofMap &,
506  const unsigned int ,
507  const Elem* )
508  { }
509 
510  template <>
512  DofMap &,
513  const unsigned int ,
514  const Elem* )
515  { }
516 #endif // LIBMESH_ENABLE_AMR
517 
518 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo