libMesh
fe_subdivision_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 // Local includes
20 #include "libmesh/fe.h"
21 #include "libmesh/libmesh_logging.h"
22 #include "libmesh/fe_type.h"
23 #include "libmesh/quadrature.h"
24 #include "libmesh/face_tri3_subdivision.h"
25 #include "libmesh/fe_macro.h"
26 #include "libmesh/dense_matrix.h"
27 #include "libmesh/utility.h"
28 
29 
30 namespace libMesh
31 {
32 
34  FE<2,SUBDIVISION>(fet)
35 {
36  // Only 2D meshes in 3D space are supported
37  libmesh_assert_equal_to(LIBMESH_DIM, 3);
38 }
39 
40 
41 
43  unsigned int valence)
44 {
45  A.resize(valence + 12, valence + 12);
46 
47  // A = (S11 0; S21 S22), see Cirak et al.,
48  // Int. J. Numer. Meth. Engng. 2000; 47:2039-2072, Appendix A.2.
49 
50  // First, set the static S21 part
51  A(valence+ 1,0 ) = 0.125;
52  A(valence+ 1,1 ) = 0.375;
53  A(valence+ 1,valence ) = 0.375;
54  A(valence+ 2,0 ) = 0.0625;
55  A(valence+ 2,1 ) = 0.625;
56  A(valence+ 2,2 ) = 0.0625;
57  A(valence+ 2,valence ) = 0.0625;
58  A(valence+ 3,0 ) = 0.125;
59  A(valence+ 3,1 ) = 0.375;
60  A(valence+ 3,2 ) = 0.375;
61  A(valence+ 4,0 ) = 0.0625;
62  A(valence+ 4,1 ) = 0.0625;
63  A(valence+ 4,valence-1) = 0.0625;
64  A(valence+ 4,valence ) = 0.625;
65  A(valence+ 5,0 ) = 0.125;
66  A(valence+ 5,valence-1) = 0.375;
67  A(valence+ 5,valence ) = 0.375;
68  A(valence+ 6,1 ) = 0.375;
69  A(valence+ 6,valence ) = 0.125;
70  A(valence+ 7,1 ) = 0.375;
71  A(valence+ 8,1 ) = 0.375;
72  A(valence+ 8,2 ) = 0.125;
73  A(valence+ 9,1 ) = 0.125;
74  A(valence+ 9,valence ) = 0.375;
75  A(valence+10,valence ) = 0.375;
76  A(valence+11,valence-1) = 0.125;
77  A(valence+11,valence ) = 0.375;
78 
79  // Next, set the static S22 part
80  A(valence+ 1,valence+1) = 0.125;
81  A(valence+ 2,valence+1) = 0.0625;
82  A(valence+ 2,valence+2) = 0.0625;
83  A(valence+ 2,valence+3) = 0.0625;
84  A(valence+ 3,valence+3) = 0.125;
85  A(valence+ 4,valence+1) = 0.0625;
86  A(valence+ 4,valence+4) = 0.0625;
87  A(valence+ 4,valence+5) = 0.0625;
88  A(valence+ 5,valence+5) = 0.125;
89  A(valence+ 6,valence+1) = 0.375;
90  A(valence+ 6,valence+2) = 0.125;
91  A(valence+ 7,valence+1) = 0.125;
92  A(valence+ 7,valence+2) = 0.375;
93  A(valence+ 7,valence+3) = 0.125;
94  A(valence+ 8,valence+2) = 0.125;
95  A(valence+ 8,valence+3) = 0.375;
96  A(valence+ 9,valence+1) = 0.375;
97  A(valence+ 9,valence+4) = 0.125;
98  A(valence+10,valence+1) = 0.125;
99  A(valence+10,valence+4) = 0.375;
100  A(valence+10,valence+5) = 0.125;
101  A(valence+11,valence+4) = 0.125;
102  A(valence+11,valence+5) = 0.375;
103 
104  // Last, set the S11 part: first row
105  std::vector<Real> weights;
106  loop_subdivision_mask(weights, valence);
107  for (unsigned int i = 0; i <= valence; ++i)
108  A(0,i) = weights[i];
109 
110  // second row
111  A(1,0) = 0.375;
112  A(1,1) = 0.375;
113  A(1,2) = 0.125;
114  A(1,valence) = 0.125;
115 
116  // third to second-to-last rows
117  for (unsigned int i = 2; i < valence; ++i)
118  {
119  A(i,0 ) = 0.375;
120  A(i,i-1) = 0.125;
121  A(i,i ) = 0.375;
122  A(i,i+1) = 0.125;
123  }
124 
125  // last row
126  A(valence,0) = 0.375;
127  A(valence,1) = 0.125;
128  A(valence,valence-1) = 0.125;
129  A(valence,valence ) = 0.375;
130 }
131 
132 
133 
134 Real FESubdivision::regular_shape(const unsigned int i,
135  const Real v,
136  const Real w)
137 {
138  // These are the 12 quartic box splines, see Cirak et al.,
139  // Int. J. Numer. Meth. Engng. 2000; 47:2039-2072, Appendix A.1.
140 
141  const Real u = 1 - v - w;
142  libmesh_assert_less_equal(0, v);
143  libmesh_assert_less_equal(0, w);
144  libmesh_assert_less_equal(0, u);
145 
146  using Utility::pow;
147  const Real factor = 1. / 12;
148 
149  switch (i)
150  {
151  case 0:
152  return factor*(pow<4>(u) + 2*u*u*u*v);
153  case 1:
154  return factor*(pow<4>(u) + 2*u*u*u*w);
155  case 2:
156  return factor*(pow<4>(u) + 2*u*u*u*w + 6*u*u*u*v + 6*u*u*v*w + 12*u*u*v*v + 6*u*v*v*w + 6*u*v*v*v +
157  2*v*v*v*w + pow<4>(v));
158  case 3:
159  return factor*(6*pow<4>(u) + 24*u*u*u*w + 24*u*u*w*w + 8*u*w*w*w + pow<4>(w) + 24*u*u*u*v +
160  60*u*u*v*w + 36*u*v*w*w + 6*v*w*w*w + 24*u*u*v*v + 36*u*v*v*w + 12*v*v*w*w + 8*u*v*v*v +
161  6*v*v*v*w + pow<4>(v));
162  case 4:
163  return factor*(pow<4>(u) + 6*u*u*u*w + 12*u*u*w*w + 6*u*w*w*w + pow<4>(w) + 2*u*u*u*v + 6*u*u*v*w +
164  6*u*v*w*w + 2*v*w*w*w);
165  case 5:
166  return factor*(2*u*v*v*v + pow<4>(v));
167  case 6:
168  return factor*(pow<4>(u) + 6*u*u*u*w + 12*u*u*w*w + 6*u*w*w*w + pow<4>(w) + 8*u*u*u*v + 36*u*u*v*w +
169  36*u*v*w*w + 8*v*w*w*w + 24*u*u*v*v + 60*u*v*v*w + 24*v*v*w*w + 24*u*v*v*v + 24*v*v*v*w + 6*pow<4>(v));
170  case 7:
171  return factor*(pow<4>(u) + 8*u*u*u*w + 24*u*u*w*w + 24*u*w*w*w + 6*pow<4>(w) + 6*u*u*u*v + 36*u*u*v*w +
172  60*u*v*w*w + 24*v*w*w*w + 12*u*u*v*v + 36*u*v*v*w + 24*v*v*w*w + 6*u*v*v*v + 8*v*v*v*w + pow<4>(v));
173  case 8:
174  return factor*(2*u*w*w*w + pow<4>(w));
175  case 9:
176  return factor*(2*v*v*v*w + pow<4>(v));
177  case 10:
178  return factor*(2*u*w*w*w + pow<4>(w) + 6*u*v*w*w + 6*v*w*w*w + 6*u*v*v*w + 12*v*v*w*w + 2*u*v*v*v +
179  6*v*v*v*w + pow<4>(v));
180  case 11:
181  return factor*(pow<4>(w) + 2*v*w*w*w);
182 
183  default:
184  libmesh_error_msg("Invalid i = " << i);
185  }
186 
187  libmesh_error_msg("We'll never get here!");
188  return 0.;
189 }
190 
191 
192 
194  const unsigned int j,
195  const Real v,
196  const Real w)
197 {
198  const Real u = 1 - v - w;
199  const Real factor = 1. / 12;
200 
201  switch (j) // j=0: xi-directional derivative, j=1: eta-directional derivative
202  {
203  case 0: // xi derivatives
204  {
205  switch (i) // shape function number
206  {
207  case 0:
208  return factor*(-6*v*u*u - 2*u*u*u);
209  case 1:
210  return factor*(-4*u*u*u - 6*u*u*w);
211  case 2:
212  return factor*(-2*v*v*v - 6*v*v*u + 6*v*u*u + 2*u*u*u);
213  case 3:
214  return factor*(-4*v*v*v - 24*v*v*u - 24*v*u*u - 18*v*v*w - 48*v*u*w - 12*u*u*w -
215  12*v*w*w - 12*u*w*w - 2*w*w*w);
216  case 4:
217  return factor*(-6*v*u*u - 2*u*u*u - 12*v*u*w-12*u*u*w - 6*v*w*w - 18*u*w*w - 4*w*w*w);
218  case 5:
219  return factor*(2*v*v*v + 6*v*v*u);
220  case 6:
221  return factor*(24*v*v*u + 24*v*u*u + 4*u*u*u + 12*v*v*w + 48*v*u*w + 18*u*u*w +
222  12*v*w*w + 12*u*w*w + 2*w*w*w);
223  case 7:
224  return factor*(-2*v*v*v - 6*v*v*u + 6*v*u*u + 2*u*u*u - 12*v*v*w + 12*u*u*w -
225  12*v*w*w + 12*u*w*w);
226  case 8:
227  return -w*w*w/6;
228  case 9:
229  return factor*(4*v*v*v + 6*v*v*w);
230  case 10:
231  return factor*(2*v*v*v + 6*v*v*u + 12*v*v*w + 12*v*u*w + 18*v*w*w + 6*u*w*w + 4*w*w*w);
232  case 11:
233  return w*w*w/6;
234  default:
235  libmesh_error_msg("Invalid i = " << i);
236  }
237  }
238  case 1: // eta derivatives
239  {
240  switch (i) // shape function number
241  {
242  case 0:
243  return factor*(-6*v*u*u - 4*u*u*u);
244  case 1:
245  return factor*(-2*u*u*u - 6*u*u*w);
246  case 2:
247  return factor*(-4*v*v*v - 18*v*v*u - 12*v*u*u - 2*u*u*u - 6*v*v*w - 12*v*u*w -
248  6*u*u*w);
249  case 3:
250  return factor*(-2*v*v*v-12*v*v*u - 12*v*u*u - 12*v*v*w - 48*v*u*w - 24*u*u*w -
251  18*v*w*w - 24*u*w*w - 4*w*w*w);
252  case 4:
253  return factor*(2*u*u*u + 6*u*u*w - 6*u*w*w - 2*w*w*w);
254  case 5:
255  return -v*v*v/6;
256  case 6:
257  return factor*(12*v*v*u + 12*v*u*u + 2*u*u*u - 12*v*v*w + 6*u*u*w - 12*v*w*w -
258  6*u*w*w - 2*w*w*w);
259  case 7:
260  return factor*(2*v*v*v + 12*v*v*u + 18*v*u*u + 4*u*u*u + 12*v*v*w + 48*v*u*w +
261  24*u*u*w + 12*v*w*w + 24*u*w*w);
262  case 8:
263  return factor*(6*u*w*w + 2*w*w*w);
264  case 9:
265  return v*v*v/6;
266  case 10:
267  return factor*(4*v*v*v + 6*v*v*u + 18*v*v*w + 12*v*u*w + 12*v*w*w + 6*u*w*w +
268  2*w*w*w);
269  case 11:
270  return factor*(6*v*w*w + 4*w*w*w);
271  default:
272  libmesh_error_msg("Invalid i = " << i);
273  }
274  }
275  default:
276  libmesh_error_msg("Invalid j = " << j);
277  }
278 
279  libmesh_error_msg("We'll never get here!");
280  return 0.;
281 }
282 
283 
284 
286  const unsigned int j,
287  const Real v,
288  const Real w)
289 {
290  const Real u = 1 - v - w;
291  const Real factor = 1. / 12;
292 
293  switch (j)
294  {
295  case 0: // xi-xi derivative
296  {
297  switch (i) // shape function number
298  {
299  case 0:
300  return v*u;
301  case 1:
302  return u*u + u*w;
303  case 2:
304  return -2*v*u;
305  case 3:
306  return v*v - 2*u*u + v*w - 2*u*w;
307  case 4:
308  return v*u + v*w + u*w + w*w;
309  case 5:
310  return v*u;
311  case 6:
312  return factor*(-24*v*v + 12*u*u - 24*v*w + 12*u*w);
313  case 7:
314  return -2*v*u - 2*v*w - 2*u*w - 2*w*w;
315  case 8:
316  return 0.;
317  case 9:
318  return v*v + v*w;
319  case 10:
320  return v*u + v*w + u*w + w*w;
321  case 11:
322  return 0.;
323  default:
324  libmesh_error_msg("Invalid i = " << i);
325  }
326  }
327  case 1: //eta-xi derivative
328  {
329  switch (i)
330  {
331  case 0:
332  return factor*(12*v*u + 6*u*u);
333  case 1:
334  return factor*(6*u*u + 12*u*w);
335  case 2:
336  return factor*(6*v*v - 12*v*u - 6*u*u);
337  case 3:
338  return factor*(6*v*v - 12*u*u + 24*v*w + 6*w*w);
339  case 4:
340  return factor*(-6*u*u - 12*u*w + 6*w*w);
341  case 5:
342  return -v*v/2.;
343  case 6:
344  return factor*(-12*v*v + 6*u*u - 24*v*w - 12*u*w - 6*w*w);
345  case 7:
346  return factor*(-6*v*v - 12*v*u + 6*u*u - 24*v*w - 12*w*w);
347  case 8:
348  return -w*w/2.;
349  case 9:
350  return v*v/2.;
351  case 10:
352  return factor*(6*v*v + 12*v*u + 24*v*w + 12*u*w + 6*w*w);
353  case 11:
354  return w*w/2.;
355  default:
356  libmesh_error_msg("Invalid i = " << i);
357  }
358  }
359  case 2: // eta-eta derivative
360  {
361  switch (i)
362  {
363  case 0:
364  return v*u + u*u;
365  case 1:
366  return u*w;
367  case 2:
368  return v*v + v*u + v*w + u*w;
369  case 3:
370  return -2*v*u - 2*u*u + v*w + w*w;
371  case 4:
372  return -2*u*w;
373  case 5:
374  return 0.;
375  case 6:
376  return -2*v*v - 2*v*u - 2*v*w - 2*u*w;
377  case 7:
378  return v*u + u*u - 2*v*w - 2*w*w;
379  case 8:
380  return u*w;
381  case 9:
382  return 0.;
383  case 10:
384  return v*v + v*u + v*w + u*w;
385  case 11:
386  return v*w + w*w;
387  default:
388  libmesh_error_msg("Invalid i = " << i);
389  }
390  }
391  default:
392  libmesh_error_msg("Invalid j = " << j);
393  }
394 
395  libmesh_error_msg("We'll never get here!");
396  return 0.;
397 }
398 
399 
400 
401 void FESubdivision::loop_subdivision_mask(std::vector<Real> & weights,
402  const unsigned int valence)
403 {
404  libmesh_assert_greater(valence, 0);
405  const Real cs = std::cos(2 * libMesh::pi / valence);
406  const Real nb_weight = (0.625 - Utility::pow<2>(0.375 + 0.25 * cs)) / valence;
407  weights.resize(1 + valence, nb_weight);
408  weights[0] = 1.0 - valence * nb_weight;
409 }
410 
411 
412 
413 void FESubdivision::init_shape_functions(const std::vector<Point> & qp,
414  const Elem * elem)
415 {
416  libmesh_assert(elem);
417  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
418  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
419 
420  LOG_SCOPE("init_shape_functions()", "FESubdivision");
421 
422  calculations_started = true;
423 
424  const unsigned int valence = sd_elem->get_ordered_valence(0);
425  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
426  const unsigned int n_approx_shape_functions = valence + 6;
427 
428  // resize the vectors to hold current data
429  phi.resize (n_approx_shape_functions);
430  dphi.resize (n_approx_shape_functions);
431  dphidxi.resize (n_approx_shape_functions);
432  dphideta.resize (n_approx_shape_functions);
433 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
434  d2phi.resize (n_approx_shape_functions);
435  d2phidxi2.resize (n_approx_shape_functions);
436  d2phidxideta.resize(n_approx_shape_functions);
437  d2phideta2.resize (n_approx_shape_functions);
438 #endif
439 
440  for (unsigned int i = 0; i < n_approx_shape_functions; ++i)
441  {
442  phi[i].resize (n_qp);
443  dphi[i].resize (n_qp);
444  dphidxi[i].resize (n_qp);
445  dphideta[i].resize (n_qp);
446 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
447  d2phi[i].resize (n_qp);
448  d2phidxi2[i].resize (n_qp);
449  d2phidxideta[i].resize(n_qp);
450  d2phideta2[i].resize (n_qp);
451 #endif
452  }
453 
454  // Renumbering of the shape functions
455  static const unsigned int cvi[12] = {3,6,2,0,1,4,7,10,9,5,11,8};
456 
457  if (valence == 6) // This means that all vertices are regular, i.e. we have 12 shape functions
458  {
459  for (unsigned int i = 0; i < n_approx_shape_functions; ++i)
460  {
461  for (unsigned int p = 0; p < n_qp; ++p)
462  {
463  phi[i][p] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, cvi[i], qp[p]);
464  dphidxi[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 0, qp[p]);
465  dphideta[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 1, qp[p]);
466  dphi[i][p](0) = dphidxi[i][p];
467  dphi[i][p](1) = dphideta[i][p];
468 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
469  d2phidxi2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 0, qp[p]);
470  d2phidxideta[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 1, qp[p]);
471  d2phideta2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 2, qp[p]);
472  d2phi[i][p](0,0) = d2phidxi2[i][p];
473  d2phi[i][p](0,1) = d2phi[i][p](1,0) = d2phidxideta[i][p];
474  d2phi[i][p](1,1) = d2phideta2[i][p];
475 #endif
476  }
477  }
478  }
479  else // vertex 0 is irregular by construction of the mesh
480  {
481  static const Real eps = 1e-10;
482 
483  // temporary values
484  std::vector<Real> tphi(12);
485  std::vector<Real> tdphidxi(12);
486  std::vector<Real> tdphideta(12);
487 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
488  std::vector<Real> td2phidxi2(12);
489  std::vector<Real> td2phidxideta(12);
490  std::vector<Real> td2phideta2(12);
491 #endif
492 
493  for (unsigned int p = 0; p < n_qp; ++p)
494  {
495  // evaluate the number of the required subdivisions
496  Real v = qp[p](0);
497  Real w = qp[p](1);
498  Real u = 1 - v - w;
499  Real min = 0, max = 0.5;
500  int n = 0;
501  while (!(u > min-eps && u < max+eps))
502  {
503  ++n;
504  min = max;
505  max += std::pow((Real)(2), -n-1);
506  }
507 
508  // transform u, v and w according to the number of subdivisions required.
509  const Real pow2 = std::pow((Real)(2), n);
510  v *= pow2;
511  w *= pow2;
512  u = 1 - v - w;
513  libmesh_assert_less(u, 0.5 + eps);
514  libmesh_assert_greater(u, -eps);
515 
516  // find out in which subdivided patch we are and setup the "selection matrix" P and the transformation Jacobian
517  // (see Int. J. Numer. Meth. Engng. 2000; 47:2039-2072, Appendix A.2.)
518  const int k = n+1;
519  Real jfac; // the additional factor per derivative order
520  DenseMatrix<Real> P(12, valence+12);
521  if (v > 0.5 - eps)
522  {
523  v = 2*v - 1;
524  w = 2*w;
525  jfac = std::pow((Real)(2), k);
526  P( 0,2 ) = 1;
527  P( 1,0 ) = 1;
528  P( 2,valence+3) = 1;
529  P( 3,1 ) = 1;
530  P( 4,valence ) = 1;
531  P( 5,valence+8) = 1;
532  P( 6,valence+2) = 1;
533  P( 7,valence+1) = 1;
534  P( 8,valence+4) = 1;
535  P( 9,valence+7) = 1;
536  P(10,valence+6) = 1;
537  P(11,valence+9) = 1;
538  }
539  else if (w > 0.5 - eps)
540  {
541  v = 2*v;
542  w = 2*w - 1;
543  jfac = std::pow((Real)(2), k);
544  P( 0,0 ) = 1;
545  P( 1,valence- 1) = 1;
546  P( 2,1 ) = 1;
547  P( 3,valence ) = 1;
548  P( 4,valence+ 5) = 1;
549  P( 5,valence+ 2) = 1;
550  P( 6,valence+ 1) = 1;
551  P( 7,valence+ 4) = 1;
552  P( 8,valence+11) = 1;
553  P( 9,valence+ 6) = 1;
554  P(10,valence+ 9) = 1;
555  P(11,valence+10) = 1;
556  }
557  else
558  {
559  v = 1 - 2*v;
560  w = 1 - 2*w;
561  jfac = std::pow((Real)(-2), k);
562  P( 0,valence+9) = 1;
563  P( 1,valence+6) = 1;
564  P( 2,valence+4) = 1;
565  P( 3,valence+1) = 1;
566  P( 4,valence+2) = 1;
567  P( 5,valence+5) = 1;
568  P( 6,valence ) = 1;
569  P( 7,1 ) = 1;
570  P( 8,valence+3) = 1;
571  P( 9,valence-1) = 1;
572  P(10,0 ) = 1;
573  P(11,2 ) = 1;
574  }
575 
576  u = 1 - v - w;
577  if ((u > 1 + eps) || (u < -eps))
578  libmesh_error_msg("SUBDIVISION irregular patch: u is outside valid range!");
579 
581  init_subdivision_matrix(A, valence);
582 
583  // compute P*A^k
584  if (k > 1)
585  {
586  DenseMatrix<Real> Acopy(A);
587  for (int e = 1; e < k; ++e)
588  A.right_multiply(Acopy);
589  }
590  P.right_multiply(A);
591 
592  const Point transformed_p(v,w);
593 
594  for (unsigned int i = 0; i < 12; ++i)
595  {
596  tphi[i] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, i, transformed_p);
597  tdphidxi[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 0, transformed_p);
598  tdphideta[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 1, transformed_p);
599 
600 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
601  td2phidxi2[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 0, transformed_p);
602  td2phidxideta[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 1, transformed_p);
603  td2phideta2[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 2, transformed_p);
604 #endif
605  }
606 
607  // Finally, we can compute the irregular shape functions as the product of P
608  // and the regular shape functions:
609  Real sum1, sum2, sum3, sum4, sum5, sum6;
610  for (unsigned int j = 0; j < n_approx_shape_functions; ++j)
611  {
612  sum1 = sum2 = sum3 = sum4 = sum5 = sum6 = 0;
613  for (unsigned int i = 0; i < 12; ++i)
614  {
615  sum1 += P(i,j) * tphi[i];
616  sum2 += P(i,j) * tdphidxi[i];
617  sum3 += P(i,j) * tdphideta[i];
618 
619 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
620  sum4 += P(i,j) * td2phidxi2[i];
621  sum5 += P(i,j) * td2phidxideta[i];
622  sum6 += P(i,j) * td2phideta2[i];
623 #endif
624  }
625  phi[j][p] = sum1;
626  dphidxi[j][p] = sum2 * jfac;
627  dphideta[j][p] = sum3 * jfac;
628  dphi[j][p](0) = dphidxi[j][p];
629  dphi[j][p](1) = dphideta[j][p];
630 
631 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
632  d2phidxi2[j][p] = sum4 * jfac * jfac;
633  d2phidxideta[j][p] = sum5 * jfac * jfac;
634  d2phideta2[j][p] = sum6 * jfac * jfac;
635  d2phi[j][p](0,0) = d2phidxi2[j][p];
636  d2phi[j][p](0,1) = d2phi[j][p](1,0) = d2phidxideta[j][p];
637  d2phi[j][p](1,1) = d2phideta2[j][p];
638 #endif
639  }
640  } // end quadrature loop
641  } // end irregular vertex
642 
643  // Let the FEMap use the same initialized shape functions
644  this->_fe_map->get_phi_map() = phi;
645  this->_fe_map->get_dphidxi_map() = dphidxi;
646  this->_fe_map->get_dphideta_map() = dphideta;
647 
648 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
649  this->_fe_map->get_d2phidxi2_map() = d2phidxi2;
650  this->_fe_map->get_d2phideta2_map() = d2phideta2;
651  this->_fe_map->get_d2phidxideta_map() = d2phidxideta;
652 #endif
653 }
654 
655 
656 
658 {
659  libmesh_assert(q);
660 
661  qrule = q;
662  // make sure we don't cache results from a previous quadrature rule
664  return;
665 }
666 
667 
668 
669 void FESubdivision::reinit(const Elem * elem,
670  const std::vector<Point> * const libmesh_dbg_var(pts),
671  const std::vector<Real> * const)
672 {
673  libmesh_assert(elem);
674  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
675 #ifndef NDEBUG
676  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
677 #endif
678 
679  LOG_SCOPE("reinit()", "FESubdivision");
680 
681  libmesh_assert(!sd_elem->is_ghost());
682  libmesh_assert(sd_elem->is_subdivision_updated());
683 
684  // check if vertices 1 and 2 are regular
685  libmesh_assert_equal_to(sd_elem->get_ordered_valence(1), 6);
686  libmesh_assert_equal_to(sd_elem->get_ordered_valence(2), 6);
687 
688  // We're calculating now! Time to determine what.
689  this->determine_calculations();
690 
691  // no custom quadrature support
694  qrule->init(elem->type());
695 
696  // Initialize the shape functions
697  this->init_shape_functions(this->qrule->get_points(), elem);
698 
699  // The shape functions correspond to the qrule
700  shapes_on_quadrature = true;
701 
702  // Compute the map for this element.
703  this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
704 }
705 
706 
707 
708 template <>
710  const Order order,
711  const unsigned int i,
712  const Point & p)
713 {
714  switch (order)
715  {
716  case FOURTH:
717  {
718  switch (type)
719  {
720  case TRI3SUBDIVISION:
721  libmesh_assert_less(i, 12);
722  return FESubdivision::regular_shape(i,p(0),p(1));
723  default:
724  libmesh_error_msg("ERROR: Unsupported element type!");
725  }
726  }
727  default:
728  libmesh_error_msg("ERROR: Unsupported polynomial order!");
729  }
730 
731  libmesh_error_msg("We'll never get here!");
732  return 0.;
733 }
734 
735 
736 
737 template <>
739  const Order order,
740  const unsigned int i,
741  const Point & p)
742 {
743  libmesh_assert(elem);
744  return FE<2,SUBDIVISION>::shape(elem->type(), order, i, p);
745 }
746 
747 
748 
749 template <>
751  const Order order,
752  const unsigned int i,
753  const unsigned int j,
754  const Point & p)
755 {
756  switch (order)
757  {
758  case FOURTH:
759  {
760  switch (type)
761  {
762  case TRI3SUBDIVISION:
763  libmesh_assert_less(i, 12);
764  return FESubdivision::regular_shape_deriv(i,j,p(0),p(1));
765  default:
766  libmesh_error_msg("ERROR: Unsupported element type!");
767  }
768  }
769  default:
770  libmesh_error_msg("ERROR: Unsupported polynomial order!");
771  }
772 
773  libmesh_error_msg("We'll never get here!");
774  return 0.;
775 }
776 
777 
778 
779 template <>
781  const Order order,
782  const unsigned int i,
783  const unsigned int j,
784  const Point & p)
785 {
786  libmesh_assert(elem);
787  return FE<2,SUBDIVISION>::shape_deriv(elem->type(), order, i, j, p);
788 }
789 
790 
791 
792 template <>
794  const Order order,
795  const unsigned int i,
796  const unsigned int j,
797  const Point & p)
798 {
799  switch (order)
800  {
801  case FOURTH:
802  {
803  switch (type)
804  {
805  case TRI3SUBDIVISION:
806  libmesh_assert_less(i, 12);
807  return FESubdivision::regular_shape_second_deriv(i,j,p(0),p(1));
808  default:
809  libmesh_error_msg("ERROR: Unsupported element type!");
810  }
811  }
812  default:
813  libmesh_error_msg("ERROR: Unsupported polynomial order!");
814  }
815 
816  libmesh_error_msg("We'll never get here!");
817  return 0.;
818 }
819 
820 
821 
822 template <>
824  const Order order,
825  const unsigned int i,
826  const unsigned int j,
827  const Point & p)
828 {
829  libmesh_assert(elem);
830  return FE<2,SUBDIVISION>::shape_second_deriv(elem->type(), order, i, j, p);
831 }
832 
833 
834 
835 template <>
837  const Order,
838  const std::vector<Number> & elem_soln,
839  std::vector<Number> & nodal_soln)
840 {
841  libmesh_assert(elem);
842  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
843  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
844 
845  nodal_soln.resize(3); // three nodes per element
846 
847  // Ghost nodes are auxiliary.
848  if (sd_elem->is_ghost())
849  {
850  nodal_soln[0] = 0;
851  nodal_soln[1] = 0;
852  nodal_soln[2] = 0;
853  return;
854  }
855 
856  // First node (node 0 in the element patch):
857  unsigned int j = sd_elem->local_node_number(sd_elem->get_ordered_node(0)->id());
858  nodal_soln[j] = elem_soln[0];
859 
860  // Second node (node 1 in the element patch):
861  j = sd_elem->local_node_number(sd_elem->get_ordered_node(1)->id());
862  nodal_soln[j] = elem_soln[1];
863 
864  // Third node (node 'valence' in the element patch):
865  j = sd_elem->local_node_number(sd_elem->get_ordered_node(2)->id());
866  nodal_soln[j] = elem_soln[sd_elem->get_ordered_valence(0)];
867 }
868 
869 
870 
871 // the empty template specializations below are needed to avoid
872 // linker reference errors, but should never get called
873 template <>
875  const Elem *,
876  const unsigned int,
877  const std::vector<Point> &,
878  std::vector<Point> &)
879 {
880  libmesh_not_implemented();
881 }
882 
883 template <>
885  unsigned int,
886  Real,
887  const std::vector<Point> * const,
888  const std::vector<Real> * const)
889 {
890  libmesh_not_implemented();
891 }
892 
893 template <>
895  const Point &,
896  const Real,
897  const bool)
898 {
899  libmesh_not_implemented();
900 }
901 
902 template <>
904  const std::vector<Point> &,
905  std::vector<Point> &,
906  Real,
907  bool)
908 {
909  libmesh_not_implemented();
910 }
911 
912 
913 
914 // The number of dofs per subdivision element depends on the valence of its nodes, so it is non-static
915 template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const ElemType, const Order) { libmesh_not_implemented(); return 0; }
916 
917 // Loop subdivision elements have only a single dof per node
918 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 1; }
919 
920 // Subdivision FEMs have dofs only at the nodes
921 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
922 
923 // Subdivision FEMs have dofs only at the nodes
924 template <> void FE<2,SUBDIVISION>::dofs_on_side(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di) { di.resize(0); }
925 template <> void FE<2,SUBDIVISION>::dofs_on_edge(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di) { di.resize(0); }
926 
927 // Subdivision FEMs are C^1 continuous
928 template <> FEContinuity FE<2,SUBDIVISION>::get_continuity() const { return C_ONE; }
929 
930 // Subdivision FEMs are not hierarchic
931 template <> bool FE<2,SUBDIVISION>::is_hierarchic() const { return false; }
932 
933 // Subdivision FEM shapes need reinit
934 template <> bool FE<2,SUBDIVISION>::shapes_need_reinit() const { return true; }
935 
936 } // namespace libMesh
bool calculate_d2phi
Should we calculate shape function hessians?
Definition: fe_abstract.h:544
std::vector< std::vector< OutputTensor > > d2phi
Shape function second derivative values.
Definition: fe_base.h:552
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:178
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Initializes the data structures for a quadrature rule for an element of type type.
Definition: quadrature.C:28
std::vector< std::vector< OutputShape > > dphidxi
Shape function derivatives in the xi direction.
Definition: fe_base.h:519
FESubdivision(const FEType &fet)
Constructor.
static unsigned int n_dofs(const ElemType t, const Order o)
static Real regular_shape(const unsigned int i, const Real v, const Real w)
virtual void edge_reinit(const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
Reinitializes all the physical element-dependent data based on the edge.
Definition: fe_boundary.C:260
virtual void side_map(const Elem *elem, const Elem *side, const unsigned int s, const std::vector< Point > &reference_side_points, std::vector< Point > &reference_points) libmesh_override
Computes the reference space quadrature points on the side of an element based on the side quadrature...
Definition: fe_boundary.C:353
bool calculations_started
Have calculations with this object already been started? Then all get_* functions should already have...
Definition: fe_abstract.h:529
bool shapes_on_quadrature
A flag indicating if current data structures correspond to quadrature rule points.
Definition: fe_abstract.h:590
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual bool shapes_need_reinit() const libmesh_override
virtual ElemType type() const =0
static void init_subdivision_matrix(DenseMatrix< Real > &A, unsigned int valence)
Builds the subdivision matrix A for the Loop scheme.
static void loop_subdivision_mask(std::vector< Real > &weights, const unsigned int valence)
Fills the vector weights with the weight coefficients of the Loop subdivision mask for evaluating the...
std::vector< std::vector< OutputShape > > d2phidxideta
Shape function second derivatives in the xi-eta direction.
Definition: fe_base.h:562
ElemType
Defines an enum for geometric element types.
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
static unsigned int n_dofs_at_node(const ElemType t, const Order o, const unsigned int n)
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static void dofs_on_side(const Elem *const elem, const Order o, unsigned int s, std::vector< unsigned int > &di)
Fills the vector di with the local degree of freedom indices associated with side s of element elem...
Definition: fe.C:63
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:197
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
This is at the core of this class.
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
virtual void init_shape_functions(const std::vector< Point > &qp, const Elem *elem) libmesh_override
Update the various member data fields phi, dphidxi, dphideta, dphidzeta, etc.
libmesh_assert(j)
A specific instantiation of the FEBase class.
Definition: fe.h:89
T pow(const T &x)
Definition: utility.h:195
std::vector< std::vector< OutputShape > > phi
Shape function values.
Definition: fe_base.h:499
std::vector< std::vector< OutputShape > > d2phideta2
Shape function second derivatives in the eta direction.
Definition: fe_base.h:572
virtual void right_multiply(const DenseMatrixBase< T > &M2) libmesh_override
Performs the operation: (*this) <- (*this) * M3.
Definition: dense_matrix.C:210
const unsigned int dim
The dimensionality of the object.
Definition: fe_abstract.h:523
virtual FEContinuity get_continuity() const libmesh_override
static Real regular_shape_second_deriv(const unsigned int i, const unsigned int j, const Real v, const Real w)
virtual void attach_quadrature_rule(QBase *q) libmesh_override
Provides the class with the quadrature rule, which provides the locations (on a reference element) wh...
virtual bool is_hierarchic() const libmesh_override
static Real regular_shape_deriv(const unsigned int i, const unsigned int j, const Real v, const Real w)
QBase * qrule
A pointer to the quadrature rule employed.
Definition: fe_abstract.h:584
PetscErrorCode Vec Mat libmesh_dbg_var(j)
static unsigned int n_dofs_per_elem(const ElemType t, const Order o)
void determine_calculations()
Determine which values are to be calculated, for both the FE itself and for the FEMap.
double pow(double a, int b)
std::vector< std::vector< OutputGradient > > dphi
Shape function derivative values.
Definition: fe_base.h:504
FEContinuity
defines an enum for finite element types to libmesh_assert a certain level (or type? Hcurl?) of continuity.
const std::vector< Real > & get_weights() const
Definition: quadrature.h:140
UniquePtr< FEMap > _fe_map
Definition: fe_abstract.h:517
static void dofs_on_edge(const Elem *const elem, const Order o, unsigned int e, std::vector< unsigned int > &di)
Fills the vector di with the local degree of freedom indices associated with edge e of element elem...
Definition: fe.C:89
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
static PetscErrorCode Mat * A
static void nodal_soln(const Elem *elem, const Order o, const std::vector< Number > &elem_soln, std::vector< Number > &nodal_soln)
Build the nodal soln from the element soln.
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
const std::vector< Point > & get_points() const
Definition: quadrature.h:128
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
std::vector< std::vector< OutputShape > > d2phidxi2
Shape function second derivatives in the xi direction.
Definition: fe_base.h:557
FEType fe_type
The finite element type for this object.
Definition: fe_abstract.h:567
ElemType elem_type
The element type the current data structures are set up for.
Definition: fe_abstract.h:573
long double min(long double a, double b)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
The QBase class provides the basic functionality from which various quadrature rules can be derived...
Definition: quadrature.h:53
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
static Point inverse_map(const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_map.C:1566
const Real pi
.
Definition: libmesh.h:172
std::vector< std::vector< OutputShape > > dphideta
Shape function derivatives in the eta direction.
Definition: fe_base.h:524