Line data Source code
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/cell_prism6.h"
21 : #include "libmesh/edge_edge2.h"
22 : #include "libmesh/face_quad4.h"
23 : #include "libmesh/face_tri3.h"
24 : #include "libmesh/enum_io_package.h"
25 : #include "libmesh/enum_order.h"
26 : #include "libmesh/fe_lagrange_shape_1D.h"
27 :
28 : // C++ includes
29 : #include <array>
30 :
31 : // Helper functions for computing volume() and true_centroid()
32 : namespace {
33 :
34 : using namespace libMesh;
35 :
36 : // Templated numerical quadrature implementation function used by both
37 : // Prism6::volume() and Prism6::true_centroid() routines.
38 : template <typename T>
39 : void cell_integral(const Elem * elem, T & t);
40 :
41 : /**
42 : * Object passed to cell_integral() routine for computing the cell volume.
43 : */
44 : struct VolumeIntegral
45 : {
46 : // API
47 2544 : void accumulate_qp(Real /*xi*/, Real /*eta*/, Real /*zeta*/, Real JxW)
48 : {
49 36840 : vol += JxW;
50 2544 : };
51 :
52 : Real vol = 0.;
53 : };
54 :
55 : /**
56 : * Object passed to cell_integral() routine for computing the nodal
57 : * volumes (one value per node)
58 : */
59 : struct NodalVolumeIntegral
60 : {
61 : // API
62 19840 : void accumulate_qp(Real xi, Real eta, Real zeta, Real JxW)
63 : {
64 : // Copied from fe_lagrange_shape_3D.C
65 : static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
66 : static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};
67 :
68 : // Copied from fe_lagrange_shape_2D.C
69 119040 : auto tri3_phi = [](int i, Real x, Real y)
70 : {
71 119040 : switch(i)
72 : {
73 39680 : case 0:
74 39680 : return 1 - x - y;
75 :
76 2304 : case 1:
77 2304 : return x;
78 :
79 39680 : case 2:
80 37376 : return y;
81 :
82 0 : default:
83 0 : libmesh_error_msg("Invalid shape function index i = " << i);
84 : }
85 : };
86 :
87 138880 : for (int i=0; i<6; ++i)
88 : {
89 : Real phi =
90 119040 : tri3_phi(i1[i], xi, eta) * fe_lagrange_1D_linear_shape(i0[i], zeta);
91 :
92 119040 : V[i] += phi * JxW;
93 : }
94 19840 : };
95 :
96 : // nodal volumes
97 : std::array<Real, 6> V {};
98 : };
99 :
100 : }
101 :
102 : namespace libMesh
103 : {
104 :
105 :
106 :
107 : // ------------------------------------------------------------
108 : // Prism6 class static member initializations
109 : const int Prism6::num_nodes;
110 : const int Prism6::nodes_per_side;
111 : const int Prism6::nodes_per_edge;
112 :
113 : const unsigned int Prism6::side_nodes_map[Prism6::num_sides][Prism6::nodes_per_side] =
114 : {
115 : {0, 2, 1, 99}, // Side 0
116 : {0, 1, 4, 3}, // Side 1
117 : {1, 2, 5, 4}, // Side 2
118 : {2, 0, 3, 5}, // Side 3
119 : {3, 4, 5, 99} // Side 4
120 : };
121 :
122 : const unsigned int Prism6::side_elems_map[Prism6::num_sides][Prism6::nodes_per_side] =
123 : {
124 : {0, 1, 2, 3}, // Side 0
125 : {0, 1, 4, 5}, // Side 1
126 : {1, 2, 5, 6}, // Side 2
127 : {0, 2, 4, 6}, // Side 3
128 : {4, 5, 6, 7} // Side 4
129 : };
130 :
131 : const unsigned int Prism6::edge_nodes_map[Prism6::num_edges][Prism6::nodes_per_edge] =
132 : {
133 : {0, 1}, // Edge 0
134 : {1, 2}, // Edge 1
135 : {0, 2}, // Edge 2
136 : {0, 3}, // Edge 3
137 : {1, 4}, // Edge 4
138 : {2, 5}, // Edge 5
139 : {3, 4}, // Edge 6
140 : {4, 5}, // Edge 7
141 : {3, 5} // Edge 8
142 : };
143 :
144 : // ------------------------------------------------------------
145 : // Prism6 class member functions
146 :
147 3589426 : bool Prism6::is_vertex(const unsigned int libmesh_dbg_var(n)) const
148 : {
149 154908 : libmesh_assert_not_equal_to (n, invalid_uint);
150 3589426 : return true;
151 : }
152 :
153 0 : bool Prism6::is_edge(const unsigned int) const
154 : {
155 0 : return false;
156 : }
157 :
158 :
159 0 : bool Prism6::is_face(const unsigned int) const
160 : {
161 0 : return false;
162 : }
163 :
164 129810 : bool Prism6::is_node_on_side(const unsigned int n,
165 : const unsigned int s) const
166 : {
167 10700 : libmesh_assert_less (s, n_sides());
168 10700 : return std::find(std::begin(side_nodes_map[s]),
169 129810 : std::end(side_nodes_map[s]),
170 129810 : n) != std::end(side_nodes_map[s]);
171 : }
172 :
173 : std::vector<unsigned>
174 40211072 : Prism6::nodes_on_side(const unsigned int s) const
175 : {
176 3586750 : libmesh_assert_less(s, n_sides());
177 40211072 : auto trim = (s > 0 && s < 4) ? 0 : 1;
178 40211072 : return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
179 : }
180 :
181 : std::vector<unsigned>
182 9279 : Prism6::nodes_on_edge(const unsigned int e) const
183 : {
184 738 : libmesh_assert_less(e, n_edges());
185 10017 : return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
186 : }
187 :
188 1759482 : bool Prism6::is_node_on_edge(const unsigned int n,
189 : const unsigned int e) const
190 : {
191 65196 : libmesh_assert_less (e, n_edges());
192 65196 : return std::find(std::begin(edge_nodes_map[e]),
193 1759482 : std::end(edge_nodes_map[e]),
194 1759482 : n) != std::end(edge_nodes_map[e]);
195 : }
196 :
197 :
198 :
199 548898 : bool Prism6::has_affine_map() const
200 : {
201 : // Make sure z edges are affine
202 107624 : Point v = this->point(3) - this->point(0);
203 645898 : if (!v.relative_fuzzy_equals(this->point(4) - this->point(1), affine_tol) ||
204 584886 : !v.relative_fuzzy_equals(this->point(5) - this->point(2), affine_tol))
205 66456 : return false;
206 48368 : return true;
207 : }
208 :
209 :
210 :
211 12131119 : Order Prism6::default_order() const
212 : {
213 12131119 : return FIRST;
214 : }
215 :
216 :
217 :
218 175186 : std::unique_ptr<Elem> Prism6::build_side_ptr (const unsigned int i)
219 : {
220 13965 : libmesh_assert_less (i, this->n_sides());
221 :
222 175186 : std::unique_ptr<Elem> face;
223 :
224 175186 : switch (i)
225 : {
226 68394 : case 0: // the triangular face at z=-1
227 : case 4: // the triangular face at z=1
228 : {
229 68394 : face = std::make_unique<Tri3>();
230 68394 : break;
231 : }
232 106792 : case 1: // the quad face at y=0
233 : case 2: // the other quad face
234 : case 3: // the quad face at x=0
235 : {
236 106792 : face = std::make_unique<Quad4>();
237 106792 : break;
238 : }
239 0 : default:
240 0 : libmesh_error_msg("Invalid side i = " << i);
241 : }
242 :
243 : // Set the nodes
244 807536 : for (auto n : face->node_index_range())
245 682762 : face->set_node(n, this->node_ptr(Prism6::side_nodes_map[i][n]));
246 :
247 175186 : face->set_interior_parent(this);
248 161221 : face->inherit_data_from(*this);
249 :
250 175186 : return face;
251 0 : }
252 :
253 :
254 :
255 7835 : void Prism6::build_side_ptr (std::unique_ptr<Elem> & side,
256 : const unsigned int i)
257 : {
258 7835 : this->side_ptr(side, i);
259 7835 : side->set_interior_parent(this);
260 6993 : side->inherit_data_from(*this);
261 7835 : }
262 :
263 :
264 :
265 1611247 : std::unique_ptr<Elem> Prism6::build_edge_ptr (const unsigned int i)
266 : {
267 1611247 : return this->simple_build_edge_ptr<Edge2,Prism6>(i);
268 : }
269 :
270 :
271 :
272 0 : void Prism6::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
273 : {
274 0 : this->simple_build_edge_ptr<Prism6>(edge, i, EDGE2);
275 0 : }
276 :
277 :
278 :
279 23040 : void Prism6::connectivity(const unsigned int libmesh_dbg_var(sc),
280 : const IOPackage iop,
281 : std::vector<dof_id_type> & conn) const
282 : {
283 1920 : libmesh_assert(_nodes);
284 1920 : libmesh_assert_less (sc, this->n_sub_elem());
285 1920 : libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
286 :
287 23040 : switch (iop)
288 : {
289 23040 : case TECPLOT:
290 : {
291 23040 : conn.resize(8);
292 24960 : conn[0] = this->node_id(0)+1;
293 23040 : conn[1] = this->node_id(1)+1;
294 23040 : conn[2] = this->node_id(2)+1;
295 23040 : conn[3] = this->node_id(2)+1;
296 23040 : conn[4] = this->node_id(3)+1;
297 23040 : conn[5] = this->node_id(4)+1;
298 23040 : conn[6] = this->node_id(5)+1;
299 23040 : conn[7] = this->node_id(5)+1;
300 23040 : return;
301 : }
302 :
303 0 : case VTK:
304 : {
305 0 : conn.resize(6);
306 0 : conn[0] = this->node_id(0);
307 0 : conn[1] = this->node_id(2);
308 0 : conn[2] = this->node_id(1);
309 0 : conn[3] = this->node_id(3);
310 0 : conn[4] = this->node_id(5);
311 0 : conn[5] = this->node_id(4);
312 0 : return;
313 : }
314 :
315 0 : default:
316 0 : libmesh_error_msg("Unsupported IO package " << iop);
317 : }
318 : }
319 :
320 :
321 :
322 : #ifdef LIBMESH_ENABLE_AMR
323 :
324 : const Real Prism6::_embedding_matrix[Prism6::num_children][Prism6::num_nodes][Prism6::num_nodes] =
325 : {
326 : // embedding matrix for child 0
327 : {
328 : // 0 1 2 3 4 5
329 : { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
330 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 1
331 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
332 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 3
333 : { .25, .25, 0.0, .25, .25, 0.0}, // 4
334 : { .25, 0.0, .25, .25, 0.0, .25} // 5
335 : },
336 :
337 : // embedding matrix for child 1
338 : {
339 : // 0 1 2 3 4 5
340 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
341 : { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
342 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 2
343 : { .25, .25, 0.0, .25, .25, 0.0}, // 3
344 : { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 4
345 : { 0.0, .25, .25, 0.0, .25, .25} // 5
346 : },
347 :
348 : // embedding matrix for child 2
349 : {
350 : // 0 1 2 3 4 5
351 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 0
352 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
353 : { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
354 : { .25, 0.0, .25, .25, 0.0, .25}, // 3
355 : { 0.0, .25, .25, 0.0, .25, .25}, // 4
356 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5} // 5
357 : },
358 :
359 : // embedding matrix for child 3
360 : {
361 : // 0 1 2 3 4 5
362 : { 0.5, 0.5, 0.0, 0.0, 0.0, 0.0}, // 0
363 : { 0.0, 0.5, 0.5, 0.0, 0.0, 0.0}, // 1
364 : { 0.5, 0.0, 0.5, 0.0, 0.0, 0.0}, // 2
365 : { .25, .25, 0.0, .25, .25, 0.0}, // 3
366 : { 0.0, .25, .25, 0.0, .25, .25}, // 4
367 : { .25, 0.0, .25, .25, 0.0, .25} // 5
368 : },
369 :
370 : // embedding matrix for child 4
371 : {
372 : // 0 1 2 3 4 5
373 : { 0.5, 0.0, 0.0, 0.5, 0.0, 0.0}, // 0
374 : { .25, .25, 0.0, .25, .25, 0.0}, // 1
375 : { .25, 0.0, .25, .25, 0.0, .25}, // 2
376 : { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 3
377 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 4
378 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
379 : },
380 :
381 : // embedding matrix for child 5
382 : {
383 : // 0 1 2 3 4 5
384 : { .25, .25, 0.0, .25, .25, 0.0}, // 0
385 : { 0.0, 0.5, 0.0, 0.0, 0.5, 0.0}, // 1
386 : { 0.0, .25, .25, 0.0, .25, .25}, // 2
387 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
388 : { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 4
389 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5} // 5
390 : },
391 :
392 : // embedding matrix for child 6
393 : {
394 : // 0 1 2 3 4 5
395 : { .25, 0.0, .25, .25, 0.0, .25}, // 0
396 : { 0.0, .25, .25, 0.0, .25, .25}, // 1
397 : { 0.0, 0.0, 0.5, 0.0, 0.0, 0.5}, // 2
398 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5}, // 3
399 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
400 : { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0} // 5
401 : },
402 :
403 : // embedding matrix for child 7
404 : {
405 : // 0 1 2 3 4 5
406 : { .25, .25, 0.0, .25, .25, 0.0}, // 0
407 : { 0.0, .25, .25, 0.0, .25, .25}, // 1
408 : { .25, 0.0, .25, .25, 0.0, .25}, // 2
409 : { 0.0, 0.0, 0.0, 0.5, 0.5, 0.0}, // 3
410 : { 0.0, 0.0, 0.0, 0.0, 0.5, 0.5}, // 4
411 : { 0.0, 0.0, 0.0, 0.5, 0.0, 0.5} // 5
412 : }
413 : };
414 :
415 : #endif
416 :
417 :
418 :
419 2480 : Point Prism6::true_centroid () const
420 : {
421 2480 : NodalVolumeIntegral nvi;
422 2480 : cell_integral(this, nvi);
423 :
424 : // Compute centroid
425 144 : Point cp;
426 144 : Real vol = 0.;
427 :
428 17360 : for (int i=0; i<6; ++i)
429 : {
430 15744 : cp += this->point(i) * nvi.V[i];
431 14880 : vol += nvi.V[i];
432 : }
433 :
434 2480 : return cp / vol;
435 : }
436 :
437 4605 : Real Prism6::volume () const
438 : {
439 : // This specialization is good for Lagrange mappings only in general
440 4605 : if (this->mapping_type() != LAGRANGE_MAP)
441 0 : return this->Elem::volume();
442 :
443 4605 : VolumeIntegral vi;
444 4605 : cell_integral(this, vi);
445 4605 : return vi.vol;
446 : }
447 :
448 : BoundingBox
449 230869 : Prism6::loose_bounding_box () const
450 : {
451 230869 : return Elem::loose_bounding_box();
452 : }
453 :
454 :
455 : void
456 4716 : Prism6::permute(unsigned int perm_num)
457 : {
458 420 : libmesh_assert_less (perm_num, 6);
459 4716 : const unsigned int side = perm_num % 2;
460 4716 : const unsigned int rotate = perm_num / 2;
461 :
462 11214 : for (unsigned int i = 0; i != rotate; ++i)
463 : {
464 6498 : swap3nodes(0,1,2);
465 5916 : swap3nodes(3,4,5);
466 5916 : swap3neighbors(1,2,3);
467 : }
468 :
469 4716 : switch (side) {
470 48 : case 0:
471 48 : break;
472 4140 : case 1:
473 4140 : swap2nodes(1,3);
474 4140 : swap2nodes(0,4);
475 4140 : swap2nodes(2,5);
476 372 : swap2neighbors(0,4);
477 372 : swap2neighbors(2,3);
478 372 : break;
479 0 : default:
480 0 : libmesh_error();
481 : }
482 :
483 4716 : }
484 :
485 :
486 : void
487 576 : Prism6::flip(BoundaryInfo * boundary_info)
488 : {
489 48 : libmesh_assert(boundary_info);
490 :
491 576 : swap2nodes(0,1);
492 576 : swap2nodes(3,4);
493 48 : swap2neighbors(2,3);
494 576 : swap2boundarysides(2,3,boundary_info);
495 576 : swap2boundaryedges(0,1,boundary_info);
496 576 : swap2boundaryedges(3,4,boundary_info);
497 576 : swap2boundaryedges(7,8,boundary_info);
498 576 : }
499 :
500 :
501 : ElemType
502 40159392 : Prism6::side_type (const unsigned int s) const
503 : {
504 3583098 : libmesh_assert_less (s, 5);
505 40159392 : if (s == 0 || s == 4)
506 2976 : return TRI3;
507 3582846 : return QUAD4;
508 : }
509 :
510 : Point
511 355 : Prism6::side_vertex_average_normal(const unsigned int s) const
512 : {
513 10 : libmesh_assert_less (s, 5);
514 10 : libmesh_assert_equal_to(this->mapping_type(), LAGRANGE_MAP);
515 :
516 355 : switch (s)
517 : {
518 142 : case 0: // the triangular face at z=-1
519 : case 4: // the triangular face at z=1
520 : {
521 150 : const Point n1 = this->point(side_nodes_map[s][1]) -
522 142 : this->point(side_nodes_map[s][0]);
523 150 : const Point n2 = this->point(side_nodes_map[s][2]) -
524 4 : this->point(side_nodes_map[s][1]);
525 138 : const Point pointing_out = n1.cross(n2);
526 142 : return pointing_out.unit();
527 : }
528 6 : case 1: // the quad face at y=0
529 : case 2: // the other quad face
530 : case 3: // the quad face at x=0
531 : {
532 : // At the side vertex average, things simplify a bit
533 : // We get the side "plane" normal at all vertices, then average them
534 6 : Point normal;
535 225 : Point current_edge = this->point(side_nodes_map[s][1]) -
536 213 : this->point(side_nodes_map[s][0]);
537 1065 : for (auto i : make_range(4))
538 : {
539 900 : const Point next_edge = this->point(side_nodes_map[s][(i + 2) % 4]) -
540 852 : this->point(side_nodes_map[s][(i + 1) % 4]);
541 828 : const Point normal_at_vertex = current_edge.cross(next_edge);
542 24 : normal += normal_at_vertex;
543 828 : current_edge = next_edge;
544 : }
545 6 : normal /= 4;
546 213 : return normal.unit();
547 : }
548 0 : default:
549 0 : libmesh_error_msg("Invalid side s = " << s);
550 : }
551 : }
552 :
553 : } // namespace libMesh
554 :
555 : namespace // anonymous
556 : {
557 :
558 : template <typename T>
559 7085 : void cell_integral(const Elem * elem, T & t)
560 : {
561 : // Conveient references to our points.
562 : const Point
563 920 : &x0 = elem->point(0), &x1 = elem->point(1), &x2 = elem->point(2),
564 462 : &x3 = elem->point(3), &x4 = elem->point(4), &x5 = elem->point(5);
565 :
566 : // constant and zeta terms only. These are copied directly from a
567 : // Python script.
568 924 : Point dx_dxi[2] =
569 : {
570 462 : -x0/2 + x1/2 - x3/2 + x4/2, // constant
571 462 : x0/2 - x1/2 - x3/2 + x4/2, // zeta
572 : };
573 :
574 : // constant and zeta terms only. These are copied directly from a
575 : // Python script.
576 924 : Point dx_deta[2] =
577 : {
578 462 : -x0/2 + x2/2 - x3/2 + x5/2, // constant
579 462 : x0/2 - x2/2 - x3/2 + x5/2, // zeta
580 : };
581 :
582 : // Constant, xi, and eta terms
583 1386 : Point dx_dzeta[3] =
584 : {
585 462 : -x0/2 + x3/2, // constant
586 462 : x0/2 - x2/2 - x3/2 + x5/2, // eta
587 462 : x0/2 - x1/2 - x3/2 + x4/2 // xi
588 : };
589 :
590 : // The quadrature rule for the Prism6 is a tensor product between a
591 : // four-point TRI3 rule (in xi, eta) and a two-point EDGE2 rule (in
592 : // zeta) which is capable of integrating cubics exactly.
593 :
594 : // Number of points in the 2D quadrature rule.
595 462 : const int N2D = 4;
596 :
597 : static const Real w2D[N2D] =
598 : {
599 : 1.5902069087198858469718450103758e-01_R,
600 : 9.0979309128011415302815498962418e-02_R,
601 : 1.5902069087198858469718450103758e-01_R,
602 : 9.0979309128011415302815498962418e-02_R
603 : };
604 :
605 : static const Real xi[N2D] =
606 : {
607 : 1.5505102572168219018027159252941e-01_R,
608 : 6.4494897427831780981972840747059e-01_R,
609 : 1.5505102572168219018027159252941e-01_R,
610 : 6.4494897427831780981972840747059e-01_R
611 : };
612 :
613 : static const Real eta[N2D] =
614 : {
615 : 1.7855872826361642311703513337422e-01_R,
616 : 7.5031110222608118177475598324603e-02_R,
617 : 6.6639024601470138670269327409637e-01_R,
618 : 2.8001991549907407200279599420481e-01_R
619 : };
620 :
621 : // Number of points in the 1D quadrature rule. The weights of the
622 : // 1D quadrature rule are equal to 1.
623 462 : const int N1D = 2;
624 :
625 : // Points of the 1D quadrature rule
626 : static const Real zeta[N1D] =
627 : {
628 : -std::sqrt(Real(3))/3,
629 : std::sqrt(Real(3))/3.
630 : };
631 :
632 35425 : for (int i=0; i<N2D; ++i)
633 : {
634 : // dx_dzeta depends only on the 2D quadrature rule points.
635 5512 : Point dx_dzeta_q = dx_dzeta[0] + eta[i]*dx_dzeta[1] + xi[i]*dx_dzeta[2];
636 :
637 85020 : for (int j=0; j<N1D; ++j)
638 : {
639 : // dx_dxi and dx_deta only depend on the 1D quadrature rule points.
640 : Point
641 7360 : dx_dxi_q = dx_dxi[0] + zeta[j]*dx_dxi[1],
642 3696 : dx_deta_q = dx_deta[0] + zeta[j]*dx_deta[1];
643 :
644 : // Compute t-specific contributions at current qp
645 60344 : t.accumulate_qp(xi[i], eta[i], zeta[j], w2D[i] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q));
646 : }
647 : }
648 7085 : }
649 :
650 : }
|
