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