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