libMesh
Loading...
Searching...
No Matches
fe_subdivision_2D.C
Go to the documentation of this file.
1// The libMesh Finite Element Library.
2// Copyright (C) 2002-2026 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
30namespace libMesh
31{
32
33
34LIBMESH_DEFAULT_VECTORIZED_FE(2,SUBDIVISION)
35
36
37FESubdivision::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
46void 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
138Real 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
194Real 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
284Real 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
398void 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
410void 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
656void 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
669void 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
708template <>
709Real 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
734template <>
735Real 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
747template <>
748Real 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
761template <>
762Real 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
788template <>
789Real 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
802template <>
803Real 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
818template <>
819Real 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
845template <>
846Real 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
860template <>
861Real 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
876template <>
877void 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
916template <>
917void 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
926template <>
927void 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
936template <>
937Point FE<2,SUBDIVISION>::inverse_map(const Elem *,
938 const Point &,
939 const Real,
940 const bool)
941{
942 libmesh_not_implemented();
943}
944
945template <>
946void 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
958template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const ElemType, const Order) { libmesh_not_implemented(); return 0; }
959template <> 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
962template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 1; }
963template <> 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
966template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
967template <> 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
970template <> void FE<2,SUBDIVISION>::dofs_on_side(const Elem * const, const Order, unsigned int, std::vector<unsigned int> & di, bool) { di.resize(0); }
971template <> 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
974template <> FEContinuity FE<2,SUBDIVISION>::get_continuity() const { return C_ONE; }
975
976// Subdivision FEMs are not hierarchic
977template <> bool FE<2,SUBDIVISION>::is_hierarchic() const { return false; }
978
979// Subdivision FEM shapes need reinit
980template <> bool FE<2,SUBDIVISION>::shapes_need_reinit() const { return true; }
981
982LIBMESH_FE_SIDE_NODAL_SOLN_DIM(SUBDIVISION, 2)
983
984} // namespace libMesh
dim
unsigned int dim
Definition adaptivity_ex3.C:114
libMesh::DenseMatrix
Defines a dense matrix for use in Finite Element-type computations.
Definition dense_matrix.h:72
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::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::DofObject::id
dof_id_type id() const
Definition dof_object.h:819
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Elem::p_level
unsigned int p_level() const
Definition elem.h:3125
libMesh::FESubdivision::FESubdivision
FESubdivision(const FEType &fet)
Constructor.
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::FEType::order
OrderWrapper order
The approximation order of the element (at 0 p-refinement level).
Definition fe_type.h:203
libMesh::FE
A specific instantiation of the FEBase class.
Definition fe.h:128
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::QBase
The QBase class provides the basic functionality from which various quadrature rules can be derived.
Definition quadrature.h:62
libMesh::Tri3Subdivision
The Tri3Subdivision element is a three-noded subdivision surface shell element used in mechanics calc...
Definition face_tri3_subdivision.h:41
libMesh::Tri3Subdivision::get_ordered_valence
unsigned int get_ordered_valence(unsigned int node_id) const
Definition face_tri3_subdivision.C:116
libMesh::Tri3Subdivision::get_ordered_node
Node * get_ordered_node(unsigned int node_id) const
Definition face_tri3_subdivision.C:124
libMesh::Tri3Subdivision::local_node_number
unsigned int local_node_number(unsigned int node_id) const
Definition face_tri3_subdivision.C:110
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::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::pi
const Real pi
.
Definition libmesh.h:292
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
Generated on Fri Aug 21 2026 13:32:29 for libMesh by doxygen 1.9.8