libMesh
fe_subdivision_2D.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 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  elem_type = INVALID_ELEM;
663  return;
664 }
665 
666 
667 
668 void FESubdivision::reinit(const Elem * elem,
669  const std::vector<Point> * const libmesh_dbg_var(pts),
670  const std::vector<Real> * const)
671 {
672  libmesh_assert(elem);
673  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
674 #ifndef NDEBUG
675  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
676 #endif
677 
678  LOG_SCOPE("reinit()", "FESubdivision");
679 
680  libmesh_assert(!sd_elem->is_ghost());
681  libmesh_assert(sd_elem->is_subdivision_updated());
682 
683  // check if vertices 1 and 2 are regular
684  libmesh_assert_equal_to(sd_elem->get_ordered_valence(1), 6);
685  libmesh_assert_equal_to(sd_elem->get_ordered_valence(2), 6);
686 
687  // We're calculating now! Time to determine what.
688  this->determine_calculations();
689 
690  // no custom quadrature support
691  libmesh_assert(pts == nullptr);
692  libmesh_assert(qrule);
693  qrule->init(elem->type());
694 
695  // Initialize the shape functions
696  this->init_shape_functions(this->qrule->get_points(), elem);
697 
698  // The shape functions correspond to the qrule
699  shapes_on_quadrature = true;
700 
701  // Compute the map for this element.
702  this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
703 }
704 
705 
706 
707 template <>
708 Real FE<2,SUBDIVISION>::shape(const ElemType type,
709  const Order order,
710  const unsigned int i,
711  const Point & p)
712 {
713  switch (order)
714  {
715  case FOURTH:
716  {
717  switch (type)
718  {
719  case TRI3SUBDIVISION:
720  libmesh_assert_less(i, 12);
721  return FESubdivision::regular_shape(i,p(0),p(1));
722  default:
723  libmesh_error_msg("ERROR: Unsupported element type == " << Utility::enum_to_string(type));
724  }
725  }
726  default:
727  libmesh_error_msg("ERROR: Unsupported polynomial order == " << order);
728  }
729 }
730 
731 
732 
733 template <>
734 Real FE<2,SUBDIVISION>::shape(const Elem * elem,
735  const Order order,
736  const unsigned int i,
737  const Point & p,
738  const bool add_p_level)
739 {
740  libmesh_assert(elem);
741  const Order totalorder =
742  static_cast<Order>(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 =
756  static_cast<Order>(fet.order+add_p_level*elem->p_level());
757  return FE<2,SUBDIVISION>::shape(elem->type(), totalorder, i, p);
758 }
759 
760 
761 
762 template <>
763 Real FE<2,SUBDIVISION>::shape_deriv(const ElemType type,
764  const Order order,
765  const unsigned int i,
766  const unsigned int j,
767  const Point & p)
768 {
769  switch (order)
770  {
771  case FOURTH:
772  {
773  switch (type)
774  {
775  case TRI3SUBDIVISION:
776  libmesh_assert_less(i, 12);
777  return FESubdivision::regular_shape_deriv(i,j,p(0),p(1));
778  default:
779  libmesh_error_msg("ERROR: Unsupported element type == " << Utility::enum_to_string(type));
780  }
781  }
782  default:
783  libmesh_error_msg("ERROR: Unsupported polynomial order == " << order);
784  }
785 }
786 
787 
788 
789 template <>
790 Real FE<2,SUBDIVISION>::shape_deriv(const Elem * elem,
791  const Order order,
792  const unsigned int i,
793  const unsigned int j,
794  const Point & p,
795  const bool add_p_level)
796 {
797  libmesh_assert(elem);
798  const Order totalorder =
799  static_cast<Order>(order+add_p_level*elem->p_level());
800  return FE<2,SUBDIVISION>::shape_deriv(elem->type(), totalorder, i, j, p);
801 }
802 
803 
804 template <>
805 Real FE<2,SUBDIVISION>::shape_deriv(const FEType fet,
806  const Elem * elem,
807  const unsigned int i,
808  const unsigned int j,
809  const Point & p,
810  const bool add_p_level)
811 {
812  libmesh_assert(elem);
813  const Order totalorder =
814  static_cast<Order>(fet.order+add_p_level*elem->p_level());
815  return FE<2,SUBDIVISION>::shape_deriv(elem->type(), totalorder, i, j, p);
816 }
817 
818 
819 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
820 
821 template <>
822 Real FE<2,SUBDIVISION>::shape_second_deriv(const ElemType type,
823  const Order order,
824  const unsigned int i,
825  const unsigned int j,
826  const Point & p)
827 {
828  switch (order)
829  {
830  case FOURTH:
831  {
832  switch (type)
833  {
834  case TRI3SUBDIVISION:
835  libmesh_assert_less(i, 12);
836  return FESubdivision::regular_shape_second_deriv(i,j,p(0),p(1));
837  default:
838  libmesh_error_msg("ERROR: Unsupported element type == " << Utility::enum_to_string(type));
839  }
840  }
841  default:
842  libmesh_error_msg("ERROR: Unsupported polynomial order == " << order);
843  }
844 }
845 
846 
847 
848 template <>
849 Real FE<2,SUBDIVISION>::shape_second_deriv(const Elem * elem,
850  const Order order,
851  const unsigned int i,
852  const unsigned int j,
853  const Point & p,
854  const bool add_p_level)
855 {
856  libmesh_assert(elem);
857  const Order totalorder =
858  static_cast<Order>(order+add_p_level*elem->p_level());
859  return FE<2,SUBDIVISION>::shape_second_deriv(elem->type(), totalorder, i, j, p);
860 }
861 
862 
863 
864 template <>
865 Real FE<2,SUBDIVISION>::shape_second_deriv(const FEType fet,
866  const Elem * elem,
867  const unsigned int i,
868  const unsigned int j,
869  const Point & p,
870  const bool add_p_level)
871 {
872  libmesh_assert(elem);
873  const Order totalorder =
874  static_cast<Order>(fet.order+add_p_level*elem->p_level());
875  return FE<2,SUBDIVISION>::shape_second_deriv(elem->type(), totalorder, i, j, p);
876 }
877 
878 
879 #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
880 
881 template <>
882 void FE<2,SUBDIVISION>::nodal_soln(const Elem * elem,
883  const Order,
884  const std::vector<Number> & elem_soln,
885  std::vector<Number> & nodal_soln,
886  const bool /*add_p_level*/)
887 {
888  libmesh_assert(elem);
889  libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
890  const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
891 
892  nodal_soln.resize(3); // three nodes per element
893 
894  // Ghost nodes are auxiliary.
895  if (sd_elem->is_ghost())
896  {
897  nodal_soln[0] = 0;
898  nodal_soln[1] = 0;
899  nodal_soln[2] = 0;
900  return;
901  }
902 
903  // First node (node 0 in the element patch):
904  unsigned int j = sd_elem->local_node_number(sd_elem->get_ordered_node(0)->id());
905  nodal_soln[j] = elem_soln[0];
906 
907  // Second node (node 1 in the element patch):
908  j = sd_elem->local_node_number(sd_elem->get_ordered_node(1)->id());
909  nodal_soln[j] = elem_soln[1];
910 
911  // Third node (node 'valence' in the element patch):
912  j = sd_elem->local_node_number(sd_elem->get_ordered_node(2)->id());
913  nodal_soln[j] = elem_soln[sd_elem->get_ordered_valence(0)];
914 }
915 
916 
917 
918 // the empty template specializations below are needed to avoid
919 // linker reference errors, but should never get called
920 template <>
921 void FE<2,SUBDIVISION>::side_map(const Elem *,
922  const Elem *,
923  const unsigned int,
924  const std::vector<Point> &,
925  std::vector<Point> &)
926 {
927  libmesh_not_implemented();
928 }
929 
930 template <>
931 void FE<2,SUBDIVISION>::edge_reinit(Elem const *,
932  unsigned int,
933  Real,
934  const std::vector<Point> * const,
935  const std::vector<Real> * const)
936 {
937  libmesh_not_implemented();
938 }
939 
940 template <>
941 Point FE<2,SUBDIVISION>::inverse_map(const Elem *,
942  const Point &,
943  const Real,
944  const bool)
945 {
946  libmesh_not_implemented();
947 }
948 
949 template <>
950 void FE<2,SUBDIVISION>::inverse_map(const Elem *,
951  const std::vector<Point> &,
952  std::vector<Point> &,
953  Real,
954  bool)
955 {
956  libmesh_not_implemented();
957 }
958 
959 
960 
961 // The number of dofs per subdivision element depends on the valence of its nodes, so it is non-static
962 template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const ElemType, const Order) { libmesh_not_implemented(); return 0; }
963 
964 // Loop subdivision elements have only a single dof per node
965 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 1; }
966 
967 // Subdivision FEMs have dofs only at the nodes
968 template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
969 
970 // Subdivision FEMs have dofs only at the nodes
971 template <> void FE<2,SUBDIVISION>::dofs_on_side(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di, bool) { di.resize(0); }
972 template <> void FE<2,SUBDIVISION>::dofs_on_edge(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di, bool) { di.resize(0); }
973 
974 // Subdivision FEMs are C^1 continuous
975 template <> FEContinuity FE<2,SUBDIVISION>::get_continuity() const { return C_ONE; }
976 
977 // Subdivision FEMs are not hierarchic
978 template <> bool FE<2,SUBDIVISION>::is_hierarchic() const { return false; }
979 
980 // Subdivision FEM shapes need reinit
981 template <> bool FE<2,SUBDIVISION>::shapes_need_reinit() const { return true; }
982 
983 LIBMESH_FE_SIDE_NODAL_SOLN_DIM(SUBDIVISION, 2)
984 
985 } // namespace libMesh
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition: fe_type.h:182
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:2945
libMesh::FEType::order
OrderWrapper order
The approximation order of the element.
Definition: fe_type.h:201
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:143
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:53
libMesh::pi
const Real pi
.
Definition: libmesh.h:274
Generated on Thu Apr 18 2024 15:17:53 for libMesh by   doxygen 1.8.14