Line data Source code
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 0 : LIBMESH_DEFAULT_VECTORIZED_FE(2,SUBDIVISION)
35 :
36 :
37 71 : FESubdivision::FESubdivision(const FEType & fet) :
38 71 : FE<2,SUBDIVISION>(fet)
39 : {
40 : // Only 2D meshes in 3D space are supported
41 : libmesh_assert_equal_to(LIBMESH_DIM, 3);
42 71 : }
43 :
44 :
45 :
46 48 : void FESubdivision::init_subdivision_matrix(DenseMatrix<Real> & A,
47 : unsigned int valence)
48 : {
49 48 : 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 52 : A(valence+ 1,0 ) = 0.125;
56 48 : A(valence+ 1,1 ) = 0.375;
57 48 : A(valence+ 1,valence ) = 0.375;
58 52 : A(valence+ 2,0 ) = 0.0625;
59 48 : A(valence+ 2,1 ) = 0.625;
60 48 : A(valence+ 2,2 ) = 0.0625;
61 48 : A(valence+ 2,valence ) = 0.0625;
62 52 : A(valence+ 3,0 ) = 0.125;
63 48 : A(valence+ 3,1 ) = 0.375;
64 48 : A(valence+ 3,2 ) = 0.375;
65 52 : A(valence+ 4,0 ) = 0.0625;
66 48 : A(valence+ 4,1 ) = 0.0625;
67 52 : A(valence+ 4,valence-1) = 0.0625;
68 48 : A(valence+ 4,valence ) = 0.625;
69 52 : A(valence+ 5,0 ) = 0.125;
70 48 : A(valence+ 5,valence-1) = 0.375;
71 48 : A(valence+ 5,valence ) = 0.375;
72 52 : A(valence+ 6,1 ) = 0.375;
73 48 : A(valence+ 6,valence ) = 0.125;
74 48 : A(valence+ 7,1 ) = 0.375;
75 52 : A(valence+ 8,1 ) = 0.375;
76 48 : A(valence+ 8,2 ) = 0.125;
77 52 : A(valence+ 9,1 ) = 0.125;
78 48 : A(valence+ 9,valence ) = 0.375;
79 48 : A(valence+10,valence ) = 0.375;
80 52 : A(valence+11,valence-1) = 0.125;
81 48 : A(valence+11,valence ) = 0.375;
82 :
83 : // Next, set the static S22 part
84 48 : A(valence+ 1,valence+1) = 0.125;
85 48 : A(valence+ 2,valence+1) = 0.0625;
86 48 : A(valence+ 2,valence+2) = 0.0625;
87 48 : A(valence+ 2,valence+3) = 0.0625;
88 48 : A(valence+ 3,valence+3) = 0.125;
89 48 : A(valence+ 4,valence+1) = 0.0625;
90 48 : A(valence+ 4,valence+4) = 0.0625;
91 48 : A(valence+ 4,valence+5) = 0.0625;
92 48 : A(valence+ 5,valence+5) = 0.125;
93 48 : A(valence+ 6,valence+1) = 0.375;
94 48 : A(valence+ 6,valence+2) = 0.125;
95 48 : A(valence+ 7,valence+1) = 0.125;
96 48 : A(valence+ 7,valence+2) = 0.375;
97 48 : A(valence+ 7,valence+3) = 0.125;
98 48 : A(valence+ 8,valence+2) = 0.125;
99 48 : A(valence+ 8,valence+3) = 0.375;
100 48 : A(valence+ 9,valence+1) = 0.375;
101 48 : A(valence+ 9,valence+4) = 0.125;
102 48 : A(valence+10,valence+1) = 0.125;
103 48 : A(valence+10,valence+4) = 0.375;
104 48 : A(valence+10,valence+5) = 0.125;
105 48 : A(valence+11,valence+4) = 0.125;
106 48 : A(valence+11,valence+5) = 0.375;
107 :
108 : // Last, set the S11 part: first row
109 4 : std::vector<Real> weights;
110 48 : loop_subdivision_mask(weights, valence);
111 288 : for (unsigned int i = 0; i <= valence; ++i)
112 280 : A(0,i) = weights[i];
113 :
114 : // second row
115 48 : A(1,0) = 0.375;
116 48 : A(1,1) = 0.375;
117 48 : A(1,2) = 0.125;
118 48 : A(1,valence) = 0.125;
119 :
120 : // third to second-to-last rows
121 144 : for (unsigned int i = 2; i < valence; ++i)
122 : {
123 96 : A(i,0 ) = 0.375;
124 104 : A(i,i-1) = 0.125;
125 96 : A(i,i ) = 0.375;
126 104 : A(i,i+1) = 0.125;
127 : }
128 :
129 : // last row
130 48 : A(valence,0) = 0.375;
131 48 : A(valence,1) = 0.125;
132 48 : A(valence,valence-1) = 0.125;
133 48 : A(valence,valence ) = 0.375;
134 48 : }
135 :
136 :
137 :
138 36864 : 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 36864 : const Real u = 1 - v - w;
146 3072 : libmesh_assert_less_equal(0, v);
147 3072 : libmesh_assert_less_equal(0, w);
148 3072 : libmesh_assert_less_equal(0, u);
149 :
150 : using Utility::pow;
151 3072 : const Real factor = 1. / 12;
152 :
153 36864 : switch (i)
154 : {
155 512 : case 0:
156 3072 : return factor*(pow<4>(u) + 2*u*u*u*v);
157 512 : case 1:
158 3072 : return factor*(pow<4>(u) + 2*u*u*u*w);
159 512 : case 2:
160 3328 : 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 3328 : 2*v*v*v*w + pow<4>(v));
162 512 : case 3:
163 3072 : 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 4096 : 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 3328 : 6*v*v*v*w + pow<4>(v));
166 512 : case 4:
167 3072 : 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 3072 : 6*u*v*w*w + 2*v*w*w*w);
169 3072 : case 5:
170 3328 : return factor*(2*u*v*v*v + pow<4>(v));
171 512 : case 6:
172 3072 : 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 3328 : 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 512 : case 7:
175 3072 : 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 3328 : 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 3072 : case 8:
178 3328 : return factor*(2*u*w*w*w + pow<4>(w));
179 3072 : case 9:
180 3328 : return factor*(2*v*v*v*w + pow<4>(v));
181 3072 : case 10:
182 3072 : 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 3328 : 6*v*v*v*w + pow<4>(v));
184 512 : case 11:
185 3072 : return factor*(pow<4>(w) + 2*v*w*w*w);
186 :
187 0 : default:
188 0 : libmesh_error_msg("Invalid i = " << i);
189 : }
190 : }
191 :
192 :
193 :
194 73728 : Real FESubdivision::regular_shape_deriv(const unsigned int i,
195 : const unsigned int j,
196 : const Real v,
197 : const Real w)
198 : {
199 73728 : const Real u = 1 - v - w;
200 6144 : const Real factor = 1. / 12;
201 :
202 73728 : switch (j) // j=0: xi-directional derivative, j=1: eta-directional derivative
203 : {
204 36864 : case 0: // xi derivatives
205 : {
206 36864 : switch (i) // shape function number
207 : {
208 3072 : case 0:
209 3072 : return factor*(-6*v*u*u - 2*u*u*u);
210 3072 : case 1:
211 3072 : return factor*(-4*u*u*u - 6*u*u*w);
212 3072 : case 2:
213 3072 : return factor*(-2*v*v*v - 6*v*v*u + 6*v*u*u + 2*u*u*u);
214 3072 : case 3:
215 3584 : 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 3072 : 12*v*w*w - 12*u*w*w - 2*w*w*w);
217 3072 : case 4:
218 3072 : 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 3072 : case 5:
220 3072 : return factor*(2*v*v*v + 6*v*v*u);
221 3072 : case 6:
222 3584 : 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 3072 : 12*v*w*w + 12*u*w*w + 2*w*w*w);
224 3072 : case 7:
225 3584 : 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 3072 : 12*v*w*w + 12*u*w*w);
227 3072 : case 8:
228 3072 : return -w*w*w/6;
229 3072 : case 9:
230 3072 : return factor*(4*v*v*v + 6*v*v*w);
231 3072 : case 10:
232 3072 : 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 3072 : case 11:
234 3072 : return w*w*w/6;
235 0 : default:
236 0 : libmesh_error_msg("Invalid i = " << i);
237 : }
238 : }
239 36864 : case 1: // eta derivatives
240 : {
241 36864 : switch (i) // shape function number
242 : {
243 3072 : case 0:
244 3072 : return factor*(-6*v*u*u - 4*u*u*u);
245 3072 : case 1:
246 3072 : return factor*(-2*u*u*u - 6*u*u*w);
247 3072 : case 2:
248 3584 : 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 3072 : 6*u*u*w);
250 3072 : case 3:
251 3584 : 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 3072 : 18*v*w*w - 24*u*w*w - 4*w*w*w);
253 3072 : case 4:
254 3072 : return factor*(2*u*u*u + 6*u*u*w - 6*u*w*w - 2*w*w*w);
255 3072 : case 5:
256 3072 : return -v*v*v/6;
257 3072 : case 6:
258 3584 : 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 3072 : 6*u*w*w - 2*w*w*w);
260 3072 : case 7:
261 3584 : 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 3072 : 24*u*u*w + 12*v*w*w + 24*u*w*w);
263 3072 : case 8:
264 3072 : return factor*(6*u*w*w + 2*w*w*w);
265 3072 : case 9:
266 3072 : return v*v*v/6;
267 3072 : case 10:
268 3584 : 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 3072 : 2*w*w*w);
270 3072 : case 11:
271 3072 : return factor*(6*v*w*w + 4*w*w*w);
272 0 : default:
273 0 : libmesh_error_msg("Invalid i = " << i);
274 : }
275 : }
276 0 : default:
277 0 : libmesh_error_msg("Invalid j = " << j);
278 : }
279 : }
280 :
281 :
282 :
283 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
284 110592 : 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 110592 : const Real u = 1 - v - w;
290 9216 : const Real factor = 1. / 12;
291 :
292 110592 : switch (j)
293 : {
294 36864 : case 0: // xi-xi derivative
295 : {
296 36864 : switch (i) // shape function number
297 : {
298 3072 : case 0:
299 3072 : return v*u;
300 3072 : case 1:
301 3072 : return u*u + u*w;
302 3072 : case 2:
303 3072 : return -2*v*u;
304 3072 : case 3:
305 3072 : return v*v - 2*u*u + v*w - 2*u*w;
306 3072 : case 4:
307 3072 : return v*u + v*w + u*w + w*w;
308 3072 : case 5:
309 3072 : return v*u;
310 3072 : case 6:
311 3072 : return factor*(-24*v*v + 12*u*u - 24*v*w + 12*u*w);
312 3072 : case 7:
313 3072 : return -2*v*u - 2*v*w - 2*u*w - 2*w*w;
314 256 : case 8:
315 256 : return 0.;
316 3072 : case 9:
317 3072 : return v*v + v*w;
318 3072 : case 10:
319 3072 : return v*u + v*w + u*w + w*w;
320 256 : case 11:
321 256 : return 0.;
322 0 : default:
323 0 : libmesh_error_msg("Invalid i = " << i);
324 : }
325 : }
326 36864 : case 1: //eta-xi derivative
327 : {
328 36864 : switch (i)
329 : {
330 3072 : case 0:
331 3072 : return factor*(12*v*u + 6*u*u);
332 3072 : case 1:
333 3072 : return factor*(6*u*u + 12*u*w);
334 3072 : case 2:
335 3072 : return factor*(6*v*v - 12*v*u - 6*u*u);
336 3072 : case 3:
337 3072 : return factor*(6*v*v - 12*u*u + 24*v*w + 6*w*w);
338 3072 : case 4:
339 3072 : return factor*(-6*u*u - 12*u*w + 6*w*w);
340 3072 : case 5:
341 3072 : return -v*v/2.;
342 3072 : case 6:
343 3072 : return factor*(-12*v*v + 6*u*u - 24*v*w - 12*u*w - 6*w*w);
344 3072 : case 7:
345 3072 : return factor*(-6*v*v - 12*v*u + 6*u*u - 24*v*w - 12*w*w);
346 3072 : case 8:
347 3072 : return -w*w/2.;
348 3072 : case 9:
349 3072 : return v*v/2.;
350 3072 : case 10:
351 3072 : return factor*(6*v*v + 12*v*u + 24*v*w + 12*u*w + 6*w*w);
352 3072 : case 11:
353 3072 : return w*w/2.;
354 0 : default:
355 0 : libmesh_error_msg("Invalid i = " << i);
356 : }
357 : }
358 36864 : case 2: // eta-eta derivative
359 : {
360 36864 : switch (i)
361 : {
362 3072 : case 0:
363 3072 : return v*u + u*u;
364 3072 : case 1:
365 3072 : return u*w;
366 3072 : case 2:
367 3072 : return v*v + v*u + v*w + u*w;
368 3072 : case 3:
369 3072 : return -2*v*u - 2*u*u + v*w + w*w;
370 3072 : case 4:
371 3072 : return -2*u*w;
372 256 : case 5:
373 256 : return 0.;
374 3072 : case 6:
375 3072 : return -2*v*v - 2*v*u - 2*v*w - 2*u*w;
376 3072 : case 7:
377 3072 : return v*u + u*u - 2*v*w - 2*w*w;
378 3072 : case 8:
379 3072 : return u*w;
380 256 : case 9:
381 256 : return 0.;
382 3072 : case 10:
383 3072 : return v*v + v*u + v*w + u*w;
384 3072 : case 11:
385 3072 : return v*w + w*w;
386 0 : default:
387 0 : libmesh_error_msg("Invalid i = " << i);
388 : }
389 : }
390 0 : default:
391 0 : libmesh_error_msg("Invalid j = " << j);
392 : }
393 : }
394 :
395 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
396 :
397 :
398 48 : void FESubdivision::loop_subdivision_mask(std::vector<Real> & weights,
399 : const unsigned int valence)
400 : {
401 4 : libmesh_assert_greater(valence, 0);
402 48 : const Real cs = std::cos(2 * libMesh::pi / valence);
403 48 : const Real nb_weight = (0.625 - Utility::pow<2>(0.375 + 0.25 * cs)) / valence;
404 48 : weights.resize(1 + valence, nb_weight);
405 48 : weights[0] = 1.0 - valence * nb_weight;
406 48 : }
407 :
408 :
409 :
410 3072 : void FESubdivision::init_shape_functions(const std::vector<Point> & qp,
411 : const Elem * elem)
412 : {
413 256 : libmesh_assert(elem);
414 256 : libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
415 256 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
416 :
417 512 : LOG_SCOPE("init_shape_functions()", "FESubdivision");
418 :
419 3072 : calculations_started = true;
420 :
421 3072 : const unsigned int valence = sd_elem->get_ordered_valence(0);
422 512 : const unsigned int n_qp = cast_int<unsigned int>(qp.size());
423 3072 : const unsigned int n_approx_shape_functions = valence + 6;
424 :
425 : // resize the vectors to hold current data
426 3072 : phi.resize (n_approx_shape_functions);
427 3072 : dphi.resize (n_approx_shape_functions);
428 3072 : dphidxi.resize (n_approx_shape_functions);
429 3072 : dphideta.resize (n_approx_shape_functions);
430 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
431 3072 : d2phi.resize (n_approx_shape_functions);
432 3072 : d2phidxi2.resize (n_approx_shape_functions);
433 3072 : d2phidxideta.resize(n_approx_shape_functions);
434 3072 : d2phideta2.resize (n_approx_shape_functions);
435 : #endif
436 :
437 39840 : for (unsigned int i = 0; i < n_approx_shape_functions; ++i)
438 : {
439 39832 : phi[i].resize (n_qp);
440 39832 : dphi[i].resize (n_qp);
441 39832 : dphidxi[i].resize (n_qp);
442 39832 : dphideta[i].resize (n_qp);
443 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
444 39832 : d2phi[i].resize (n_qp);
445 39832 : d2phidxi2[i].resize (n_qp);
446 39832 : d2phidxideta[i].resize(n_qp);
447 39832 : 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 3072 : if (valence == 6) // This means that all vertices are regular, i.e. we have 12 shape functions
455 : {
456 39312 : for (unsigned int i = 0; i < n_approx_shape_functions; ++i)
457 : {
458 72576 : for (unsigned int p = 0; p < n_qp; ++p)
459 : {
460 39312 : phi[i][p] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, cvi[i], qp[p]);
461 39312 : dphidxi[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 0, qp[p]);
462 39312 : dphideta[i][p] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, cvi[i], 1, qp[p]);
463 45360 : dphi[i][p](0) = dphidxi[i][p];
464 36288 : dphi[i][p](1) = dphideta[i][p];
465 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
466 39312 : d2phidxi2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 0, qp[p]);
467 39312 : d2phidxideta[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 1, qp[p]);
468 39312 : d2phideta2[i][p] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, cvi[i], 2, qp[p]);
469 45360 : d2phi[i][p](0,0) = d2phidxi2[i][p];
470 39312 : d2phi[i][p](0,1) = d2phi[i][p](1,0) = d2phidxideta[i][p];
471 36288 : 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 52 : std::vector<Real> tphi(12);
482 52 : std::vector<Real> tdphidxi(12);
483 52 : std::vector<Real> tdphideta(12);
484 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
485 52 : std::vector<Real> td2phidxi2(12);
486 52 : std::vector<Real> td2phidxideta(12);
487 52 : std::vector<Real> td2phideta2(12);
488 : #endif
489 :
490 96 : for (unsigned int p = 0; p < n_qp; ++p)
491 : {
492 : // evaluate the number of the required subdivisions
493 48 : Real v = qp[p](0);
494 48 : Real w = qp[p](1);
495 48 : Real u = 1 - v - w;
496 4 : Real min = 0, max = 0.5;
497 4 : int n = 0;
498 48 : while (!(u > min-eps && u < max+eps))
499 : {
500 0 : ++n;
501 0 : min = max;
502 0 : max += std::pow((Real)(2), -n-1);
503 : }
504 :
505 : // transform u, v and w according to the number of subdivisions required.
506 4 : const Real pow2 = std::pow((Real)(2), n);
507 48 : v *= pow2;
508 48 : w *= pow2;
509 4 : u = 1 - v - w;
510 4 : libmesh_assert_less(u, 0.5 + eps);
511 4 : 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 48 : const int k = n+1;
516 : Real jfac; // the additional factor per derivative order
517 56 : DenseMatrix<Real> P(12, valence+12);
518 48 : if (v > 0.5 - eps)
519 : {
520 0 : v = 2*v - 1;
521 0 : w = 2*w;
522 0 : jfac = std::pow((Real)(2), k);
523 0 : P( 0,2 ) = 1;
524 0 : P( 1,0 ) = 1;
525 0 : P( 2,valence+3) = 1;
526 0 : P( 3,1 ) = 1;
527 0 : P( 4,valence ) = 1;
528 0 : P( 5,valence+8) = 1;
529 0 : P( 6,valence+2) = 1;
530 0 : P( 7,valence+1) = 1;
531 0 : P( 8,valence+4) = 1;
532 0 : P( 9,valence+7) = 1;
533 0 : P(10,valence+6) = 1;
534 0 : P(11,valence+9) = 1;
535 : }
536 48 : else if (w > 0.5 - eps)
537 : {
538 0 : v = 2*v;
539 0 : w = 2*w - 1;
540 0 : jfac = std::pow((Real)(2), k);
541 0 : P( 0,0 ) = 1;
542 0 : P( 1,valence- 1) = 1;
543 0 : P( 2,1 ) = 1;
544 0 : P( 3,valence ) = 1;
545 0 : P( 4,valence+ 5) = 1;
546 0 : P( 5,valence+ 2) = 1;
547 0 : P( 6,valence+ 1) = 1;
548 0 : P( 7,valence+ 4) = 1;
549 0 : P( 8,valence+11) = 1;
550 0 : P( 9,valence+ 6) = 1;
551 0 : P(10,valence+ 9) = 1;
552 0 : P(11,valence+10) = 1;
553 : }
554 : else
555 : {
556 48 : v = 1 - 2*v;
557 48 : w = 1 - 2*w;
558 4 : jfac = std::pow((Real)(-2), k);
559 52 : P( 0,valence+9) = 1;
560 48 : P( 1,valence+6) = 1;
561 48 : P( 2,valence+4) = 1;
562 48 : P( 3,valence+1) = 1;
563 48 : P( 4,valence+2) = 1;
564 52 : P( 5,valence+5) = 1;
565 48 : P( 6,valence ) = 1;
566 48 : P( 7,1 ) = 1;
567 48 : P( 8,valence+3) = 1;
568 52 : P( 9,valence-1) = 1;
569 48 : P(10,0 ) = 1;
570 48 : P(11,2 ) = 1;
571 : }
572 :
573 48 : u = 1 - v - w;
574 48 : libmesh_error_msg_if((u > 1 + eps) || (u < -eps), "SUBDIVISION irregular patch: u is outside valid range!");
575 :
576 56 : DenseMatrix<Real> A;
577 48 : init_subdivision_matrix(A, valence);
578 :
579 : // compute P*A^k
580 48 : if (k > 1)
581 : {
582 0 : DenseMatrix<Real> Acopy(A);
583 0 : for (int e = 1; e < k; ++e)
584 0 : A.right_multiply(Acopy);
585 0 : }
586 48 : P.right_multiply(A);
587 :
588 4 : const Point transformed_p(v,w);
589 :
590 624 : for (unsigned int i = 0; i < 12; ++i)
591 : {
592 576 : tphi[i] = FE<2,SUBDIVISION>::shape (elem, fe_type.order, i, transformed_p);
593 576 : tdphidxi[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 0, transformed_p);
594 576 : tdphideta[i] = FE<2,SUBDIVISION>::shape_deriv (elem, fe_type.order, i, 1, transformed_p);
595 :
596 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
597 576 : td2phidxi2[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 0, transformed_p);
598 576 : td2phidxideta[i] = FE<2,SUBDIVISION>::shape_second_deriv(elem, fe_type.order, i, 1, transformed_p);
599 576 : 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 528 : for (unsigned int j = 0; j < n_approx_shape_functions; ++j)
607 : {
608 40 : sum1 = sum2 = sum3 = sum4 = sum5 = sum6 = 0;
609 6240 : for (unsigned int i = 0; i < 12; ++i)
610 : {
611 5760 : sum1 += P(i,j) * tphi[i];
612 5760 : sum2 += P(i,j) * tdphidxi[i];
613 5760 : sum3 += P(i,j) * tdphideta[i];
614 :
615 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
616 5760 : sum4 += P(i,j) * td2phidxi2[i];
617 5760 : sum5 += P(i,j) * td2phidxideta[i];
618 6240 : sum6 += P(i,j) * td2phideta2[i];
619 : #endif
620 : }
621 520 : phi[j][p] = sum1;
622 520 : dphidxi[j][p] = sum2 * jfac;
623 520 : dphideta[j][p] = sum3 * jfac;
624 520 : dphi[j][p](0) = dphidxi[j][p];
625 480 : dphi[j][p](1) = dphideta[j][p];
626 :
627 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
628 520 : d2phidxi2[j][p] = sum4 * jfac * jfac;
629 520 : d2phidxideta[j][p] = sum5 * jfac * jfac;
630 520 : d2phideta2[j][p] = sum6 * jfac * jfac;
631 520 : d2phi[j][p](0,0) = d2phidxi2[j][p];
632 480 : d2phi[j][p](0,1) = d2phi[j][p](1,0) = d2phidxideta[j][p];
633 480 : d2phi[j][p](1,1) = d2phideta2[j][p];
634 : #endif
635 : }
636 40 : } // end quadrature loop
637 : } // end irregular vertex
638 :
639 : // Let the FEMap use the same initialized shape functions
640 3072 : this->_fe_map->get_phi_map() = phi;
641 3072 : this->_fe_map->get_dphidxi_map() = dphidxi;
642 3072 : this->_fe_map->get_dphideta_map() = dphideta;
643 :
644 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
645 3072 : this->_fe_map->get_d2phidxi2_map() = d2phidxi2;
646 3072 : this->_fe_map->get_d2phideta2_map() = d2phideta2;
647 3072 : this->_fe_map->get_d2phidxideta_map() = d2phidxideta;
648 : #endif
649 :
650 3072 : if (this->calculate_dual)
651 0 : this->init_dual_shape_functions(n_approx_shape_functions, n_qp);
652 3072 : }
653 :
654 :
655 :
656 71 : void FESubdivision::attach_quadrature_rule(QBase * q)
657 : {
658 2 : libmesh_assert(q);
659 :
660 71 : qrule = q;
661 : // make sure we don't cache results from a previous quadrature rule
662 71 : this->_elem = nullptr;
663 71 : this->_elem_type = INVALID_ELEM;
664 71 : return;
665 : }
666 :
667 :
668 :
669 3072 : void FESubdivision::reinit(const Elem * elem,
670 : const std::vector<Point> * const libmesh_dbg_var(pts),
671 : const std::vector<Real> * const)
672 : {
673 256 : libmesh_assert(elem);
674 256 : libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
675 : #ifndef NDEBUG
676 256 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
677 : #endif
678 :
679 512 : LOG_SCOPE("reinit()", "FESubdivision");
680 :
681 256 : libmesh_assert(!sd_elem->is_ghost());
682 256 : libmesh_assert(sd_elem->is_subdivision_updated());
683 :
684 : // check if vertices 1 and 2 are regular
685 256 : libmesh_assert_equal_to(sd_elem->get_ordered_valence(1), 6);
686 256 : libmesh_assert_equal_to(sd_elem->get_ordered_valence(2), 6);
687 :
688 : // We're calculating now! Time to determine what.
689 3072 : this->determine_calculations();
690 :
691 : // no custom quadrature support
692 256 : libmesh_assert(pts == nullptr);
693 256 : libmesh_assert(qrule);
694 3072 : qrule->init(*elem);
695 :
696 : // Initialize the shape functions
697 3328 : this->init_shape_functions(this->qrule->get_points(), elem);
698 :
699 : // The shape functions correspond to the qrule
700 3072 : shapes_on_quadrature = true;
701 :
702 : // Compute the map for this element.
703 3328 : this->_fe_map->compute_map (this->dim, this->qrule->get_weights(), elem, this->calculate_d2phi);
704 3072 : }
705 :
706 :
707 :
708 : template <>
709 36864 : Real FE<2,SUBDIVISION>::shape(const ElemType type,
710 : const Order order,
711 : const unsigned int i,
712 : const Point & p)
713 : {
714 36864 : switch (order)
715 : {
716 36864 : case FOURTH:
717 : {
718 36864 : switch (type)
719 : {
720 36864 : case TRI3SUBDIVISION:
721 3072 : libmesh_assert_less(i, 12);
722 36864 : return FESubdivision::regular_shape(i,p(0),p(1));
723 0 : default:
724 0 : libmesh_error_msg("ERROR: Unsupported element type == " << Utility::enum_to_string(type));
725 : }
726 : }
727 0 : default:
728 0 : libmesh_error_msg("ERROR: Unsupported polynomial order == " << order);
729 : }
730 : }
731 :
732 :
733 :
734 : template <>
735 36864 : 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 3072 : libmesh_assert(elem);
742 39936 : const Order totalorder = order + add_p_level*elem->p_level();
743 36864 : return FE<2,SUBDIVISION>::shape(elem->type(), totalorder, i, p);
744 : }
745 :
746 :
747 : template <>
748 0 : 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 0 : libmesh_assert(elem);
755 0 : const Order totalorder = fet.order + add_p_level*elem->p_level();
756 0 : return FE<2,SUBDIVISION>::shape(elem->type(), totalorder, i, p);
757 : }
758 :
759 :
760 :
761 : template <>
762 73728 : 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 73728 : switch (order)
769 : {
770 73728 : case FOURTH:
771 : {
772 73728 : switch (type)
773 : {
774 73728 : case TRI3SUBDIVISION:
775 6144 : libmesh_assert_less(i, 12);
776 73728 : return FESubdivision::regular_shape_deriv(i,j,p(0),p(1));
777 0 : default:
778 0 : libmesh_error_msg("ERROR: Unsupported element type == " << Utility::enum_to_string(type));
779 : }
780 : }
781 0 : default:
782 0 : libmesh_error_msg("ERROR: Unsupported polynomial order == " << order);
783 : }
784 : }
785 :
786 :
787 :
788 : template <>
789 73728 : 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 6144 : libmesh_assert(elem);
797 79872 : const Order totalorder = order + add_p_level*elem->p_level();
798 73728 : return FE<2,SUBDIVISION>::shape_deriv(elem->type(), totalorder, i, j, p);
799 : }
800 :
801 :
802 : template <>
803 0 : 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 0 : libmesh_assert(elem);
811 0 : const Order totalorder = fet.order + add_p_level*elem->p_level();
812 0 : 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 110592 : 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 110592 : switch (order)
826 : {
827 110592 : case FOURTH:
828 : {
829 110592 : switch (type)
830 : {
831 110592 : case TRI3SUBDIVISION:
832 9216 : libmesh_assert_less(i, 12);
833 110592 : return FESubdivision::regular_shape_second_deriv(i,j,p(0),p(1));
834 0 : default:
835 0 : libmesh_error_msg("ERROR: Unsupported element type == " << Utility::enum_to_string(type));
836 : }
837 : }
838 0 : default:
839 0 : libmesh_error_msg("ERROR: Unsupported polynomial order == " << order);
840 : }
841 : }
842 :
843 :
844 :
845 : template <>
846 110592 : 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 9216 : libmesh_assert(elem);
854 119808 : const Order totalorder = order + add_p_level*elem->p_level();
855 110592 : return FE<2,SUBDIVISION>::shape_second_deriv(elem->type(), totalorder, i, j, p);
856 : }
857 :
858 :
859 :
860 : template <>
861 0 : 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 0 : libmesh_assert(elem);
869 0 : const Order totalorder = fet.order + add_p_level*elem->p_level();
870 0 : 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 11664 : 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 972 : libmesh_assert(elem);
885 972 : libmesh_assert_equal_to(elem->type(), TRI3SUBDIVISION);
886 972 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
887 :
888 11664 : nodal_soln.resize(3); // three nodes per element
889 :
890 : // Ghost nodes are auxiliary.
891 11664 : if (sd_elem->is_ghost())
892 : {
893 2244 : nodal_soln[0] = 0;
894 2244 : nodal_soln[1] = 0;
895 2244 : nodal_soln[2] = 0;
896 2448 : return;
897 : }
898 :
899 : // First node (node 0 in the element patch):
900 9216 : unsigned int j = sd_elem->local_node_number(sd_elem->get_ordered_node(0)->id());
901 9216 : nodal_soln[j] = elem_soln[0];
902 :
903 : // Second node (node 1 in the element patch):
904 9216 : j = sd_elem->local_node_number(sd_elem->get_ordered_node(1)->id());
905 9216 : nodal_soln[j] = elem_soln[1];
906 :
907 : // Third node (node 'valence' in the element patch):
908 9216 : j = sd_elem->local_node_number(sd_elem->get_ordered_node(2)->id());
909 9216 : 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 0 : 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 0 : libmesh_not_implemented();
924 : }
925 :
926 : template <>
927 0 : 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 0 : libmesh_not_implemented();
934 : }
935 :
936 : template <>
937 0 : Point FE<2,SUBDIVISION>::inverse_map(const Elem *,
938 : const Point &,
939 : const Real,
940 : const bool)
941 : {
942 0 : libmesh_not_implemented();
943 : }
944 :
945 : template <>
946 0 : void FE<2,SUBDIVISION>::inverse_map(const Elem *,
947 : const std::vector<Point> &,
948 : std::vector<Point> &,
949 : Real,
950 : bool)
951 : {
952 0 : 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 0 : template <> unsigned int FE<2,SUBDIVISION>::n_dofs(const ElemType, const Order) { libmesh_not_implemented(); return 0; }
959 0 : 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 562734 : template <> unsigned int FE<2,SUBDIVISION>::n_dofs_at_node(const ElemType, const Order, const unsigned int) { return 1; }
963 138024 : 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 0 : template <> unsigned int FE<2,SUBDIVISION>::n_dofs_per_elem(const ElemType, const Order) { return 0; }
967 70023 : 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 0 : 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 0 : 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 0 : template <> FEContinuity FE<2,SUBDIVISION>::get_continuity() const { return C_ONE; }
975 :
976 : // Subdivision FEMs are not hierarchic
977 0 : template <> bool FE<2,SUBDIVISION>::is_hierarchic() const { return false; }
978 :
979 : // Subdivision FEM shapes need reinit
980 0 : template <> bool FE<2,SUBDIVISION>::shapes_need_reinit() const { return true; }
981 :
982 0 : LIBMESH_FE_SIDE_NODAL_SOLN_DIM(SUBDIVISION, 2)
983 :
984 : } // namespace libMesh
|
