fe_subdivision_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 
33 FESubdivision::FESubdivision(const FEType & fet) :
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 
42 void FESubdivision::init_subdivision_matrix(DenseMatrix<Real> & A,
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 
188 
189 
190 Real FESubdivision::regular_shape_deriv(const unsigned int i,
191  const unsigned int j,
192  const Real v,
193  const Real w)
194 {
195  const Real u = 1 - v - w;
196  const Real factor = 1. / 12;
197 
198  switch (j) // j=0: xi-directional derivative, j=1: eta-directional derivative
199  {
200  case 0: // xi derivatives
201  {
202  switch (i) // shape function number
203  {
204  case 0:
205  return factor*(-6*v*u*u - 2*u*u*u);
206  case 1:
207  return factor*(-4*u*u*u - 6*u*u*w);
208  case 2:
209  return factor*(-2*v*v*v - 6*v*v*u + 6*v*u*u + 2*u*u*u);
210  case 3:
211  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 -
212  12*v*w*w - 12*u*w*w - 2*w*w*w);
213  case 4:
214  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);
215  case 5:
216  return factor*(2*v*v*v + 6*v*v*u);
217  case 6:
218  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 +
219  12*v*w*w + 12*u*w*w + 2*w*w*w);
220  case 7:
221  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 -
222  12*v*w*w + 12*u*w*w);
223  case 8:
224  return -w*w*w/6;
225  case 9:
226  return factor*(4*v*v*v + 6*v*v*w);
227  case 10:
228  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);
229  case 11:
230  return w*w*w/6;
231  default:
232  libmesh_error_msg("Invalid i = " << i);
233  }
234  }
235  case 1: // eta derivatives
236  {
237  switch (i) // shape function number
238  {
239  case 0:
240  return factor*(-6*v*u*u - 4*u*u*u);
241  case 1:
242  return factor*(-2*u*u*u - 6*u*u*w);
243  case 2:
244  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 -
245  6*u*u*w);
246  case 3:
247  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 -
248  18*v*w*w - 24*u*w*w - 4*w*w*w);
249  case 4:
250  return factor*(2*u*u*u + 6*u*u*w - 6*u*w*w - 2*w*w*w);
251  case 5:
252  return -v*v*v/6;
253  case 6:
254  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 -
255  6*u*w*w - 2*w*w*w);
256  case 7:
257  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 +
258  24*u*u*w + 12*v*w*w + 24*u*w*w);
259  case 8:
260  return factor*(6*u*w*w + 2*w*w*w);
261  case 9:
262  return v*v*v/6;
263  case 10:
264  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 +
265  2*w*w*w);
266  case 11:
267  return factor*(6*v*w*w + 4*w*w*w);
268  default:
269  libmesh_error_msg("Invalid i = " << i);
270  }
271  }
272  default:
273  libmesh_error_msg("Invalid j = " << j);
274  }
275 }
276 
277 
278 
279 Real FESubdivision::regular_shape_second_deriv(const unsigned int i,
280  const unsigned int j,
281  const Real v,
282  const Real w)
283 {
284  const Real u = 1 - v - w;
285  const Real factor = 1. / 12;
286 
287  switch (j)
288  {
289  case 0: // xi-xi derivative
290  {
291  switch (i) // shape function number
292  {
293  case 0:
294  return v*u;
295  case 1:
296  return u*u + u*w;
297  case 2:
298  return -2*v*u;
299  case 3:
300  return v*v - 2*u*u + v*w - 2*u*w;
301  case 4:
302  return v*u + v*w + u*w + w*w;
303  case 5:
304  return v*u;
305  case 6:
306  return factor*(-24*v*v + 12*u*u - 24*v*w + 12*u*w);
307  case 7:
308  return -2*v*u - 2*v*w - 2*u*w - 2*w*w;
309  case 8:
310  return 0.;
311  case 9:
312  return v*v + v*w;
313  case 10:
314  return v*u + v*w + u*w + w*w;
315  case 11:
316  return 0.;
317  default:
318  libmesh_error_msg("Invalid i = " << i);
319  }
320  }
321  case 1: //eta-xi derivative
322  {
323  switch (i)
324  {
325  case 0:
326  return factor*(12*v*u + 6*u*u);
327  case 1:
328  return factor*(6*u*u + 12*u*w);
329  case 2:
330  return factor*(6*v*v - 12*v*u - 6*u*u);
331  case 3:
332  return factor*(6*v*v - 12*u*u + 24*v*w + 6*w*w);
333  case 4:
334  return factor*(-6*u*u - 12*u*w + 6*w*w);
335  case 5:
336  return -v*v/2.;
337  case 6:
338  return factor*(-12*v*v + 6*u*u - 24*v*w - 12*u*w - 6*w*w);
339  case 7:
340  return factor*(-6*v*v - 12*v*u + 6*u*u - 24*v*w - 12*w*w);
341  case 8:
342  return -w*w/2.;
343  case 9:
344  return v*v/2.;
345  case 10:
346  return factor*(6*v*v + 12*v*u + 24*v*w + 12*u*w + 6*w*w);
347  case 11:
348  return w*w/2.;
349  default:
350  libmesh_error_msg("Invalid i = " << i);
351  }
352  }
353  case 2: // eta-eta derivative
354  {
355  switch (i)
356  {
357  case 0:
358  return v*u + u*u;
359  case 1:
360  return u*w;
361  case 2:
362  return v*v + v*u + v*w + u*w;
363  case 3:
364  return -2*v*u - 2*u*u + v*w + w*w;
365  case 4:
366  return -2*u*w;
367  case 5:
368  return 0.;
369  case 6:
370  return -2*v*v - 2*v*u - 2*v*w - 2*u*w;
371  case 7:
372  return v*u + u*u - 2*v*w - 2*w*w;
373  case 8:
374  return u*w;
375  case 9:
376  return 0.;
377  case 10:
378  return v*v + v*u + v*w + u*w;
379  case 11:
380  return v*w + w*w;
381  default:
382  libmesh_error_msg("Invalid i = " << i);
383  }
384  }
385  default:
386  libmesh_error_msg("Invalid j = " << j);
387  }
388 }
389 
390 
391 
392 void FESubdivision::loop_subdivision_mask(std::vector<Real> & weights,
393  const unsigned int valence)
394 {
395  libmesh_assert_greater(valence, 0);
396  const Real cs = std::cos(2 * libMesh::pi / valence);
397  const Real nb_weight = (0.625 - Utility::pow<2>(0.375 + 0.25 * cs)) / valence;
398  weights.resize(1 + valence, nb_weight);
399  weights[0] = 1.0 - valence * nb_weight;
400 }
401 
402 
403 
404 void FESubdivision::init_shape_functions(const std::vector<Point> & qp,
405  const Elem * elem)
406 {
407  libmesh_assert(elem);
408  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
409  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
410 
411  LOG_SCOPE("init_shape_functions()", "FESubdivision");
412 
413  calculations_started = true;
414 
415  const unsigned int valence = sd_elem->get_ordered_valence(0);
416  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
417  const unsigned int n_approx_shape_functions = valence + 6;
418 
419  // resize the vectors to hold current data
420  phi.resize (n_approx_shape_functions);
421  dphi.resize (n_approx_shape_functions);
422  dphidxi.resize (n_approx_shape_functions);
423  dphideta.resize (n_approx_shape_functions);
424 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
425  d2phi.resize (n_approx_shape_functions);
426  d2phidxi2.resize (n_approx_shape_functions);
427  d2phidxideta.resize(n_approx_shape_functions);
428  d2phideta2.resize (n_approx_shape_functions);
429 #endif
430 
431  for (unsigned int i = 0; i < n_approx_shape_functions; ++i)
432  {
433  phi[i].resize (n_qp);
434  dphi[i].resize (n_qp);
435  dphidxi[i].resize (n_qp);
436  dphideta[i].resize (n_qp);
437 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
438  d2phi[i].resize (n_qp);
439  d2phidxi2[i].resize (n_qp);
440  d2phidxideta[i].resize(n_qp);
441  d2phideta2[i].resize (n_qp);
442 #endif
443  }
444 
445  // Renumbering of the shape functions
446  static const unsigned int cvi[12] = {3,6,2,0,1,4,7,10,9,5,11,8};
447 
448  if (valence == 6) // This means that all vertices are regular, i.e. we have 12 shape functions
449  {
450  for (unsigned int i = 0; i < n_approx_shape_functions; ++i)
451  {
452  for (unsigned int p = 0; p < n_qp; ++p)
453  {
454  phi[i][p] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, cvi[i], qp[p]);
455  dphidxi[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 0, qp[p]);
456  dphideta[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 1, qp[p]);
457  dphi[i][p](0) = dphidxi[i][p];
458  dphi[i][p](1) = dphideta[i][p];
459 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
460  d2phidxi2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 0, qp[p]);
461  d2phidxideta[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 1, qp[p]);
462  d2phideta2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 2, qp[p]);
463  d2phi[i][p](0,0) = d2phidxi2[i][p];
464  d2phi[i][p](0,1) = d2phi[i][p](1,0) = d2phidxideta[i][p];
465  d2phi[i][p](1,1) = d2phideta2[i][p];
466 #endif
467  }
468  }
469  }
470  else // vertex 0 is irregular by construction of the mesh
471  {
472  static const Real eps = 1e-10;
473 
474  // temporary values
475  std::vector<Real> tphi(12);
476  std::vector<Real> tdphidxi(12);
477  std::vector<Real> tdphideta(12);
478 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
479  std::vector<Real> td2phidxi2(12);
480  std::vector<Real> td2phidxideta(12);
481  std::vector<Real> td2phideta2(12);
482 #endif
483 
484  for (unsigned int p = 0; p < n_qp; ++p)
485  {
486  // evaluate the number of the required subdivisions
487  Real v = qp[p](0);
488  Real w = qp[p](1);
489  Real u = 1 - v - w;
490  Real min = 0, max = 0.5;
491  int n = 0;
492  while (!(u > min-eps && u < max+eps))
493  {
494  ++n;
495  min = max;
496  max += std::pow((Real)(2), -n-1);
497  }
498 
499  // transform u, v and w according to the number of subdivisions required.
500  const Real pow2 = std::pow((Real)(2), n);
501  v *= pow2;
502  w *= pow2;
503  u = 1 - v - w;
504  libmesh_assert_less(u, 0.5 + eps);
505  libmesh_assert_greater(u, -eps);
506 
507  // find out in which subdivided patch we are and setup the "selection matrix" P and the transformation Jacobian
508  // (see Int. J. Numer. Meth. Engng. 2000; 47:2039-2072, Appendix A.2.)
509  const int k = n+1;
510  Real jfac; // the additional factor per derivative order
511  DenseMatrix<Real> P(12, valence+12);
512  if (v > 0.5 - eps)
513  {
514  v = 2*v - 1;
515  w = 2*w;
516  jfac = std::pow((Real)(2), k);
517  P( 0,2 ) = 1;
518  P( 1,0 ) = 1;
519  P( 2,valence+3) = 1;
520  P( 3,1 ) = 1;
521  P( 4,valence ) = 1;
522  P( 5,valence+8) = 1;
523  P( 6,valence+2) = 1;
524  P( 7,valence+1) = 1;
525  P( 8,valence+4) = 1;
526  P( 9,valence+7) = 1;
527  P(10,valence+6) = 1;
528  P(11,valence+9) = 1;
529  }
530  else if (w > 0.5 - eps)
531  {
532  v = 2*v;
533  w = 2*w - 1;
534  jfac = std::pow((Real)(2), k);
535  P( 0,0 ) = 1;
536  P( 1,valence- 1) = 1;
537  P( 2,1 ) = 1;
538  P( 3,valence ) = 1;
539  P( 4,valence+ 5) = 1;
540  P( 5,valence+ 2) = 1;
541  P( 6,valence+ 1) = 1;
542  P( 7,valence+ 4) = 1;
543  P( 8,valence+11) = 1;
544  P( 9,valence+ 6) = 1;
545  P(10,valence+ 9) = 1;
546  P(11,valence+10) = 1;
547  }
548  else
549  {
550  v = 1 - 2*v;
551  w = 1 - 2*w;
552  jfac = std::pow((Real)(-2), k);
553  P( 0,valence+9) = 1;
554  P( 1,valence+6) = 1;
555  P( 2,valence+4) = 1;
556  P( 3,valence+1) = 1;
557  P( 4,valence+2) = 1;
558  P( 5,valence+5) = 1;
559  P( 6,valence ) = 1;
560  P( 7,1 ) = 1;
561  P( 8,valence+3) = 1;
562  P( 9,valence-1) = 1;
563  P(10,0 ) = 1;
564  P(11,2 ) = 1;
565  }
566 
567  u = 1 - v - w;
568  if ((u > 1 + eps) || (u < -eps))
569  libmesh_error_msg("SUBDIVISION irregular patch: u is outside valid range!");
570 
571  DenseMatrix<Real> A;
572  init_subdivision_matrix(A, valence);
573 
574  // compute P*A^k
575  if (k > 1)
576  {
577  DenseMatrix<Real> Acopy(A);
578  for (int e = 1; e < k; ++e)
579  A.right_multiply(Acopy);
580  }
581  P.right_multiply(A);
582 
583  const Point transformed_p(v,w);
584 
585  for (unsigned int i = 0; i < 12; ++i)
586  {
587  tphi[i] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, i, transformed_p);
588  tdphidxi[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 0, transformed_p);
589  tdphideta[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 1, transformed_p);
590 
591 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
592  td2phidxi2[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 0, transformed_p);
593  td2phidxideta[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 1, transformed_p);
594  td2phideta2[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 2, transformed_p);
595 #endif
596  }
597 
598  // Finally, we can compute the irregular shape functions as the product of P
599  // and the regular shape functions:
600  Real sum1, sum2, sum3, sum4, sum5, sum6;
601  for (unsigned int j = 0; j < n_approx_shape_functions; ++j)
602  {
603  sum1 = sum2 = sum3 = sum4 = sum5 = sum6 = 0;
604  for (unsigned int i = 0; i < 12; ++i)
605  {
606  sum1 += P(i,j) * tphi[i];
607  sum2 += P(i,j) * tdphidxi[i];
608  sum3 += P(i,j) * tdphideta[i];
609 
610 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
611  sum4 += P(i,j) * td2phidxi2[i];
612  sum5 += P(i,j) * td2phidxideta[i];
613  sum6 += P(i,j) * td2phideta2[i];
614 #endif
615  }
616  phi[j][p] = sum1;
617  dphidxi[j][p] = sum2 * jfac;
618  dphideta[j][p] = sum3 * jfac;
619  dphi[j][p](0) = dphidxi[j][p];
620  dphi[j][p](1) = dphideta[j][p];
621 
622 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
623  d2phidxi2[j][p] = sum4 * jfac * jfac;
624  d2phidxideta[j][p] = sum5 * jfac * jfac;
625  d2phideta2[j][p] = sum6 * jfac * jfac;
626  d2phi[j][p](0,0) = d2phidxi2[j][p];
627  d2phi[j][p](0,1) = d2phi[j][p](1,0) = d2phidxideta[j][p];
628  d2phi[j][p](1,1) = d2phideta2[j][p];
629 #endif
630  }
631  } // end quadrature loop
632  } // end irregular vertex
633 
634  // Let the FEMap use the same initialized shape functions
635  this->_fe_map->get_phi_map() = phi;
636  this->_fe_map->get_dphidxi_map() = dphidxi;
637  this->_fe_map->get_dphideta_map() = dphideta;
638 
639 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
640  this->_fe_map->get_d2phidxi2_map() = d2phidxi2;
641  this->_fe_map->get_d2phideta2_map() = d2phideta2;
642  this->_fe_map->get_d2phidxideta_map() = d2phidxideta;
643 #endif
644 }
645 
646 
647 
648 void FESubdivision::attach_quadrature_rule(QBase * q)
649 {
650  libmesh_assert(q);
651 
652  qrule = q;
653  // make sure we don't cache results from a previous quadrature rule
654  elem_type = INVALID_ELEM;
655  return;
656 }
657 
658 
659 
660 void FESubdivision::reinit(const Elem * elem,
661  const std::vector<Point> * const libmesh_dbg_var(pts),
662  const std::vector<Real> * const)
663 {
664  libmesh_assert(elem);
665  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
666 #ifndef NDEBUG
667  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
668 #endif
669 
670  LOG_SCOPE("reinit()", "FESubdivision");
671 
672  libmesh_assert(!sd_elem->is_ghost());
673  libmesh_assert(sd_elem->is_subdivision_updated());
674 
675  // check if vertices 1 and 2 are regular
676  libmesh_assert_equal_to(sd_elem->get_ordered_valence(1), 6);
677  libmesh_assert_equal_to(sd_elem->get_ordered_valence(2), 6);
678 
679  // We're calculating now! Time to determine what.
680  this->determine_calculations();
681 
682  // no custom quadrature support
683  libmesh_assert(pts == nullptr);
684  libmesh_assert(qrule);
685  qrule->init(elem->type());
686 
687  // Initialize the shape functions
688  this->init_shape_functions(this->qrule->get_points(), elem);
689 
690  // The shape functions correspond to the qrule
691  shapes_on_quadrature = true;
692 
693  // Compute the map for this element.
694  this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
695 }
696 
697 
698 
699 template <>
700 Real FE<2,SUBDIVISION>::shape(const ElemType type,
701  const Order order,
702  const unsigned int i,
703  const Point & p)
704 {
705  switch (order)
706  {
707  case FOURTH:
708  {
709  switch (type)
710  {
711  case TRI3SUBDIVISION:
712  libmesh_assert_less(i, 12);
713  return FESubdivision::regular_shape(i,p(0),p(1));
714  default:
715  libmesh_error_msg("ERROR: Unsupported element type!");
716  }
717  }
718  default:
719  libmesh_error_msg("ERROR: Unsupported polynomial order!");
720  }
721 }
722 
723 
724 
725 template <>
726 Real FE<2,SUBDIVISION>::shape(const Elem * elem,
727  const Order order,
728  const unsigned int i,
729  const Point & p)
730 {
731  libmesh_assert(elem);
732  return FE<2,SUBDIVISION>::shape(elem->type(), order, i, p);
733 }
734 
735 
736 
737 template <>
738 Real FE<2,SUBDIVISION>::shape_deriv(const ElemType type,
739  const Order order,
740  const unsigned int i,
741  const unsigned int j,
742  const Point & p)
743 {
744  switch (order)
745  {
746  case FOURTH:
747  {
748  switch (type)
749  {
750  case TRI3SUBDIVISION:
751  libmesh_assert_less(i, 12);
752  return FESubdivision::regular_shape_deriv(i,j,p(0),p(1));
753  default:
754  libmesh_error_msg("ERROR: Unsupported element type!");
755  }
756  }
757  default:
758  libmesh_error_msg("ERROR: Unsupported polynomial order!");
759  }
760 }
761 
762 
763 
764 template <>
765 Real FE<2,SUBDIVISION>::shape_deriv(const Elem * elem,
766  const Order order,
767  const unsigned int i,
768  const unsigned int j,
769  const Point & p)
770 {
771  libmesh_assert(elem);
772  return FE<2,SUBDIVISION>::shape_deriv(elem->type(), order, i, j, p);
773 }
774 
775 
776 
777 template <>
778 Real FE<2,SUBDIVISION>::shape_second_deriv(const ElemType type,
779  const Order order,
780  const unsigned int i,
781  const unsigned int j,
782  const Point & p)
783 {
784  switch (order)
785  {
786  case FOURTH:
787  {
788  switch (type)
789  {
790  case TRI3SUBDIVISION:
791  libmesh_assert_less(i, 12);
792  return FESubdivision::regular_shape_second_deriv(i,j,p(0),p(1));
793  default:
794  libmesh_error_msg("ERROR: Unsupported element type!");
795  }
796  }
797  default:
798  libmesh_error_msg("ERROR: Unsupported polynomial order!");
799  }
800 }
801 
802 
803 
804 template <>
805 Real FE<2,SUBDIVISION>::shape_second_deriv(const Elem * elem,
806  const Order order,
807  const unsigned int i,
808  const unsigned int j,
809  const Point & p)
810 {
811  libmesh_assert(elem);
812  return FE<2,SUBDIVISION>::shape_second_deriv(elem->type(), order, i, j, p);
813 }
814 
815 
816 
817 template <>
818 void FE<2,SUBDIVISION>::nodal_soln(const Elem * elem,
819  const Order,
820  const std::vector<Number> & elem_soln,
821  std::vector<Number> & nodal_soln)
822 {
823  libmesh_assert(elem);
824  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
825  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
826 
827  nodal_soln.resize(3); // three nodes per element
828 
829  // Ghost nodes are auxiliary.
830  if (sd_elem->is_ghost())
831  {
832  nodal_soln[0] = 0;
833  nodal_soln[1] = 0;
834  nodal_soln[2] = 0;
835  return;
836  }
837 
838  // First node (node 0 in the element patch):
839  unsigned int j = sd_elem->local_node_number(sd_elem->get_ordered_node(0)->id());
840  nodal_soln[j] = elem_soln[0];
841 
842  // Second node (node 1 in the element patch):
843  j = sd_elem->local_node_number(sd_elem->get_ordered_node(1)->id());
844  nodal_soln[j] = elem_soln[1];
845 
846  // Third node (node 'valence' in the element patch):
847  j = sd_elem->local_node_number(sd_elem->get_ordered_node(2)->id());
848  nodal_soln[j] = elem_soln[sd_elem->get_ordered_valence(0)];
849 }
850 
851 
852 
853 // the empty template specializations below are needed to avoid
854 // linker reference errors, but should never get called
855 template <>
856 void FE<2,SUBDIVISION>::side_map(const Elem *,
857  const Elem *,
858  const unsigned int,
859  const std::vector<Point> &,
860  std::vector<Point> &)
861 {
862  libmesh_not_implemented();
863 }
864 
865 template <>
866 void FE<2,SUBDIVISION>::edge_reinit(Elem const *,
867  unsigned int,
868  Real,
869  const std::vector<Point> * const,
870  const std::vector<Real> * const)
871 {
872  libmesh_not_implemented();
873 }
874 
875 template <>
876 Point FE<2,SUBDIVISION>::inverse_map(const Elem *,
877  const Point &,
878  const Real,
879  const bool)
880 {
881  libmesh_not_implemented();
882 }
883 
884 template <>
885 void FE<2,SUBDIVISION>::inverse_map(const Elem *,
886  const std::vector<Point> &,
887  std::vector<Point> &,
888  Real,
889  bool)
890 {
891  libmesh_not_implemented();
892 }
893 
894 
895 
896 // The number of dofs per subdivision element depends on the valence of its nodes, so it is non-static
897 template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const ElemType, const Order) { libmesh_not_implemented(); return 0; }
898 
899 // Loop subdivision elements have only a single dof per node
900 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 1; }
901 
902 // Subdivision FEMs have dofs only at the nodes
903 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
904 
905 // Subdivision FEMs have dofs only at the nodes
906 template <> void FE<2,SUBDIVISION>::dofs_on_side(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di) { di.resize(0); }
907 template <> void FE<2,SUBDIVISION>::dofs_on_edge(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di) { di.resize(0); }
908 
909 // Subdivision FEMs are C^1 continuous
910 template <> FEContinuity FE<2,SUBDIVISION>::get_continuity() const { return C_ONE; }
911 
912 // Subdivision FEMs are not hierarchic
913 template <> bool FE<2,SUBDIVISION>::is_hierarchic() const { return false; }
914 
915 // Subdivision FEM shapes need reinit
916 template <> bool FE<2,SUBDIVISION>::shapes_need_reinit() const { return true; }
917 
918 } // namespace libMesh
libMesh::FEGenericBase< FEOutputType< T >::type >::d2phi
std::vector< std::vector< OutputTensor > > d2phi
Definition: fe_base.h:551
libMesh::FEType
Manages the family, order, etc. parameters for a given FE.
Definition: fe_type.h:179
libMesh::QBase::init
virtual void init(const ElemType type=INVALID_ELEM, unsigned int p_level=0)
Definition: quadrature.C:28
libMesh::FEGenericBase< FEOutputType< T >::type >::dphidxi
std::vector< std::vector< OutputShape > > dphidxi
Definition: fe_base.h:518
libMesh::FESubdivision::FESubdivision
FESubdivision(const FEType &fet)
Definition: fe_subdivision_2D.C:33
libMesh::FESubdivision::regular_shape
static Real regular_shape(const unsigned int i, const Real v, const Real w)
Definition: fe_subdivision_2D.C:134
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::FESubdivision::init_subdivision_matrix
static void init_subdivision_matrix(DenseMatrix< Real > &A, unsigned int valence)
Definition: fe_subdivision_2D.C:42
libMesh::FESubdivision::loop_subdivision_mask
static void loop_subdivision_mask(std::vector< Real > &weights, const unsigned int valence)
Definition: fe_subdivision_2D.C:392
libMesh::FEGenericBase< FEOutputType< T >::type >::d2phidxideta
std::vector< std::vector< OutputShape > > d2phidxideta
Definition: fe_base.h:561
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::QBase::get_weights
const std::vector< Real > & get_weights() const
Definition: quadrature.h:154
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FEType::order
OrderWrapper order
Definition: fe_type.h:198
std::max
long double max(long double a, double b)
Definition: libmesh_augment_std_namespace.h:27
libMesh::FESubdivision::attach_quadrature_rule
virtual void attach_quadrature_rule(QBase *q) override
Definition: fe_subdivision_2D.C:648
dense_matrix.h
libMesh::FE
Template class which generates the different FE families and orders.
Definition: fe.h:89
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:198
libMesh::FEGenericBase< FEOutputType< T >::type >::phi
std::vector< std::vector< OutputShape > > phi
Definition: fe_base.h:498
libMesh::FEGenericBase< FEOutputType< T >::type >::d2phideta2
std::vector< std::vector< OutputShape > > d2phideta2
Definition: fe_base.h:571
libMesh::FEAbstract::dim
const unsigned int dim
Definition: fe_abstract.h:531
libMesh::FESubdivision::regular_shape_second_deriv
static Real regular_shape_second_deriv(const unsigned int i, const unsigned int j, const Real v, const Real w)
Definition: fe_subdivision_2D.C:279
libMesh::DenseMatrix::right_multiply
virtual void right_multiply(const DenseMatrixBase< T > &M2) override
Definition: dense_matrix_impl.h:213
libmesh_logging.h
libMesh::Tri3Subdivision
A surface shell element used in mechanics calculations.
Definition: face_tri3_subdivision.h:40
libMesh::FESubdivision::regular_shape_deriv
static Real regular_shape_deriv(const unsigned int i, const unsigned int j, const Real v, const Real w)
Definition: fe_subdivision_2D.C:190
libMesh::Tri3Subdivision::local_node_number
unsigned int local_node_number(unsigned int node_id) const
Definition: face_tri3_subdivision.C:106
face_tri3_subdivision.h
std::pow
double pow(double a, int b)
Definition: libmesh_augment_std_namespace.h:66
libMesh::FEGenericBase< FEOutputType< T >::type >::dphi
std::vector< std::vector< OutputGradient > > dphi
Definition: fe_base.h:503
fe_macro.h
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::QBase::get_points
const std::vector< Point > & get_points() const
Definition: quadrature.h:142
A
static PetscErrorCode Mat * A
Definition: petscdmlibmeshimpl.C:1020
libMesh::FEGenericBase< FEOutputType< T >::type >::d2phidxi2
std::vector< std::vector< OutputShape > > d2phidxi2
Definition: fe_base.h:556
fe_type.h
libMesh::FESubdivision::init_shape_functions
virtual void init_shape_functions(const std::vector< Point > &qp, const Elem *elem) override
Definition: fe_subdivision_2D.C:404
libMesh::FEAbstract::_fe_map
std::unique_ptr< FEMap > _fe_map
Definition: fe_abstract.h:525
libMesh::DenseMatrix
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54
libMesh::Elem::type
virtual ElemType type() const =0
std::min
long double min(long double a, double b)
Definition: libmesh_augment_std_namespace.h:29
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::FESubdivision::reinit
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr) override
libMesh::QBase
Base class for all quadrature families and orders.
Definition: quadrature.h:62
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::pi
const Real pi
Definition: libmesh.h:233
libMesh::FEGenericBase< FEOutputType< T >::type >::dphideta
std::vector< std::vector< OutputShape > > dphideta
Definition: fe_base.h:523
generated by doxygen