MFEM v4.7.0 Finite element discretization library
Searching...
No Matches
fe_fixed_order.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12// Fixed Order Finite Element classes
13
14#include "fe_fixed_order.hpp"
15#include "../coefficient.hpp"
16
17namespace mfem
18{
19
20using namespace std;
21
23 : NodalFiniteElement(0, Geometry::POINT, 1, 0)
24{
26 lex_ordering[0] = 0;
27 Nodes.IntPoint(0).x = 0.0;
28}
29
31 Vector &shape) const
32{
33 shape(0) = 1.;
34}
35
37 DenseMatrix &dshape) const
38{
39 // dshape is (1 x 0) - nothing to compute
40}
41
43 : NodalFiniteElement(1, Geometry::SEGMENT, 2, 1)
44{
45 Nodes.IntPoint(0).x = 0.0;
46 Nodes.IntPoint(1).x = 1.0;
47}
48
50 Vector &shape) const
51{
52 shape(0) = 1. - ip.x;
53 shape(1) = ip.x;
54}
55
57 DenseMatrix &dshape) const
58{
59 dshape(0,0) = -1.;
60 dshape(1,0) = 1.;
61}
62
64 : NodalFiniteElement(2, Geometry::TRIANGLE, 3, 1)
65{
66 Nodes.IntPoint(0).x = 0.0;
67 Nodes.IntPoint(0).y = 0.0;
68 Nodes.IntPoint(1).x = 1.0;
69 Nodes.IntPoint(1).y = 0.0;
70 Nodes.IntPoint(2).x = 0.0;
71 Nodes.IntPoint(2).y = 1.0;
72}
73
75 Vector &shape) const
76{
77 shape(0) = 1. - ip.x - ip.y;
78 shape(1) = ip.x;
79 shape(2) = ip.y;
80}
81
83 DenseMatrix &dshape) const
84{
85 dshape(0,0) = -1.; dshape(0,1) = -1.;
86 dshape(1,0) = 1.; dshape(1,1) = 0.;
87 dshape(2,0) = 0.; dshape(2,1) = 1.;
88}
89
90
92 : NodalFiniteElement(2, Geometry::SQUARE, 4, 1, FunctionSpace::Qk)
93{
94 Nodes.IntPoint(0).x = 0.0;
95 Nodes.IntPoint(0).y = 0.0;
96 Nodes.IntPoint(1).x = 1.0;
97 Nodes.IntPoint(1).y = 0.0;
98 Nodes.IntPoint(2).x = 1.0;
99 Nodes.IntPoint(2).y = 1.0;
100 Nodes.IntPoint(3).x = 0.0;
101 Nodes.IntPoint(3).y = 1.0;
102}
103
105 Vector &shape) const
106{
107 shape(0) = (1. - ip.x) * (1. - ip.y) ;
108 shape(1) = ip.x * (1. - ip.y) ;
109 shape(2) = ip.x * ip.y ;
110 shape(3) = (1. - ip.x) * ip.y ;
111}
112
114 DenseMatrix &dshape) const
115{
116 dshape(0,0) = -1. + ip.y; dshape(0,1) = -1. + ip.x ;
117 dshape(1,0) = 1. - ip.y; dshape(1,1) = -ip.x ;
118 dshape(2,0) = ip.y ; dshape(2,1) = ip.x ;
119 dshape(3,0) = -ip.y ; dshape(3,1) = 1. - ip.x ;
120}
121
123 const IntegrationPoint &ip, DenseMatrix &h) const
124{
125 h(0,0) = 0.; h(0,1) = 1.; h(0,2) = 0.;
126 h(1,0) = 0.; h(1,1) = -1.; h(1,2) = 0.;
127 h(2,0) = 0.; h(2,1) = 1.; h(2,2) = 0.;
128 h(3,0) = 0.; h(3,1) = -1.; h(3,2) = 0.;
129}
130
131
133 : NodalFiniteElement(2, Geometry::TRIANGLE, 3, 1, FunctionSpace::Pk)
134{
135 Nodes.IntPoint(0).x = 1./6.;
136 Nodes.IntPoint(0).y = 1./6.;
137 Nodes.IntPoint(1).x = 2./3.;
138 Nodes.IntPoint(1).y = 1./6.;
139 Nodes.IntPoint(2).x = 1./6.;
140 Nodes.IntPoint(2).y = 2./3.;
141}
142
144 Vector &shape) const
145{
146 const real_t x = ip.x, y = ip.y;
147
148 shape(0) = 5./3. - 2. * (x + y);
149 shape(1) = 2. * (x - 1./6.);
150 shape(2) = 2. * (y - 1./6.);
151}
152
154 DenseMatrix &dshape) const
155{
156 dshape(0,0) = -2.; dshape(0,1) = -2.;
157 dshape(1,0) = 2.; dshape(1,1) = 0.;
158 dshape(2,0) = 0.; dshape(2,1) = 2.;
159}
160
162{
163 dofs(vertex) = 2./3.;
164 dofs((vertex+1)%3) = 1./6.;
165 dofs((vertex+2)%3) = 1./6.;
166}
167
168
169// 0.5-0.5/sqrt(3) and 0.5+0.5/sqrt(3)
170const real_t GaussBiLinear2DFiniteElement::p[] =
171{ 0.2113248654051871177454256, 0.7886751345948128822545744 };
172
174 : NodalFiniteElement(2, Geometry::SQUARE, 4, 1, FunctionSpace::Qk)
175{
176 Nodes.IntPoint(0).x = p[0];
177 Nodes.IntPoint(0).y = p[0];
178 Nodes.IntPoint(1).x = p[1];
179 Nodes.IntPoint(1).y = p[0];
180 Nodes.IntPoint(2).x = p[1];
181 Nodes.IntPoint(2).y = p[1];
182 Nodes.IntPoint(3).x = p[0];
183 Nodes.IntPoint(3).y = p[1];
184}
185
187 Vector &shape) const
188{
189 const real_t x = ip.x, y = ip.y;
190
191 shape(0) = 3. * (p[1] - x) * (p[1] - y);
192 shape(1) = 3. * (x - p[0]) * (p[1] - y);
193 shape(2) = 3. * (x - p[0]) * (y - p[0]);
194 shape(3) = 3. * (p[1] - x) * (y - p[0]);
195}
196
198 DenseMatrix &dshape) const
199{
200 const real_t x = ip.x, y = ip.y;
201
202 dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
203 dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
204 dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
205 dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
206}
207
209{
210#if 1
211 dofs(vertex) = p[1]*p[1];
212 dofs((vertex+1)%4) = p[0]*p[1];
213 dofs((vertex+2)%4) = p[0]*p[0];
214 dofs((vertex+3)%4) = p[0]*p[1];
215#else
216 dofs = 1.0;
217#endif
218}
219
220
222 : NodalFiniteElement(2, Geometry::SQUARE, 3, 1, FunctionSpace::Qk)
223{
224 Nodes.IntPoint(0).x = 0.0;
225 Nodes.IntPoint(0).y = 0.0;
226 Nodes.IntPoint(1).x = 1.0;
227 Nodes.IntPoint(1).y = 0.0;
228 Nodes.IntPoint(2).x = 0.0;
229 Nodes.IntPoint(2).y = 1.0;
230}
231
233 Vector &shape) const
234{
235 shape(0) = 1. - ip.x - ip.y;
236 shape(1) = ip.x;
237 shape(2) = ip.y;
238}
239
241 DenseMatrix &dshape) const
242{
243 dshape(0,0) = -1.; dshape(0,1) = -1.;
244 dshape(1,0) = 1.; dshape(1,1) = 0.;
245 dshape(2,0) = 0.; dshape(2,1) = 1.;
246}
247
248
250 : NodalFiniteElement(1, Geometry::SEGMENT, 3, 2)
251{
252 Nodes.IntPoint(0).x = 0.0;
253 Nodes.IntPoint(1).x = 1.0;
254 Nodes.IntPoint(2).x = 0.5;
255}
256
258 Vector &shape) const
259{
260 real_t x = ip.x;
261 real_t l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
262
263 shape(0) = l1 * (-l3);
264 shape(1) = l2 * l3;
265 shape(2) = 4. * l1 * l2;
266}
267
269 DenseMatrix &dshape) const
270{
271 real_t x = ip.x;
272
273 dshape(0,0) = 4. * x - 3.;
274 dshape(1,0) = 4. * x - 1.;
275 dshape(2,0) = 4. - 8. * x;
276}
277
278
280 : NodalFiniteElement(2, Geometry::TRIANGLE, 6, 2)
281{
282 Nodes.IntPoint(0).x = 0.0;
283 Nodes.IntPoint(0).y = 0.0;
284 Nodes.IntPoint(1).x = 1.0;
285 Nodes.IntPoint(1).y = 0.0;
286 Nodes.IntPoint(2).x = 0.0;
287 Nodes.IntPoint(2).y = 1.0;
288 Nodes.IntPoint(3).x = 0.5;
289 Nodes.IntPoint(3).y = 0.0;
290 Nodes.IntPoint(4).x = 0.5;
291 Nodes.IntPoint(4).y = 0.5;
292 Nodes.IntPoint(5).x = 0.0;
293 Nodes.IntPoint(5).y = 0.5;
294}
295
297 Vector &shape) const
298{
299 real_t x = ip.x, y = ip.y;
300 real_t l1 = 1.-x-y, l2 = x, l3 = y;
301
302 shape(0) = l1 * (2. * l1 - 1.);
303 shape(1) = l2 * (2. * l2 - 1.);
304 shape(2) = l3 * (2. * l3 - 1.);
305 shape(3) = 4. * l1 * l2;
306 shape(4) = 4. * l2 * l3;
307 shape(5) = 4. * l3 * l1;
308}
309
311 DenseMatrix &dshape) const
312{
313 real_t x = ip.x, y = ip.y;
314
315 dshape(0,0) =
316 dshape(0,1) = 4. * (x + y) - 3.;
317
318 dshape(1,0) = 4. * x - 1.;
319 dshape(1,1) = 0.;
320
321 dshape(2,0) = 0.;
322 dshape(2,1) = 4. * y - 1.;
323
324 dshape(3,0) = -4. * (2. * x + y - 1.);
325 dshape(3,1) = -4. * x;
326
327 dshape(4,0) = 4. * y;
328 dshape(4,1) = 4. * x;
329
330 dshape(5,0) = -4. * y;
331 dshape(5,1) = -4. * (x + 2. * y - 1.);
332}
333
335 DenseMatrix &h) const
336{
337 h(0,0) = 4.;
338 h(0,1) = 4.;
339 h(0,2) = 4.;
340
341 h(1,0) = 4.;
342 h(1,1) = 0.;
343 h(1,2) = 0.;
344
345 h(2,0) = 0.;
346 h(2,1) = 0.;
347 h(2,2) = 4.;
348
349 h(3,0) = -8.;
350 h(3,1) = -4.;
351 h(3,2) = 0.;
352
353 h(4,0) = 0.;
354 h(4,1) = 4.;
355 h(4,2) = 0.;
356
357 h(5,0) = 0.;
358 h(5,1) = -4.;
359 h(5,2) = -8.;
360}
361
362void Quad2DFiniteElement::ProjectDelta(int vertex, Vector &dofs) const
363{
364#if 0
365 dofs = 1.;
366#else
367 dofs = 0.;
368 dofs(vertex) = 1.;
369 switch (vertex)
370 {
371 case 0: dofs(3) = 0.25; dofs(5) = 0.25; break;
372 case 1: dofs(3) = 0.25; dofs(4) = 0.25; break;
373 case 2: dofs(4) = 0.25; dofs(5) = 0.25; break;
374 }
375#endif
376}
377
378
380{ 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
381
383 : NodalFiniteElement(2, Geometry::TRIANGLE, 6, 2), A(6), D(6,2), pol(6)
384{
385 Nodes.IntPoint(0).x = p[0];
386 Nodes.IntPoint(0).y = p[0];
387 Nodes.IntPoint(1).x = 1. - 2. * p[0];
388 Nodes.IntPoint(1).y = p[0];
389 Nodes.IntPoint(2).x = p[0];
390 Nodes.IntPoint(2).y = 1. - 2. * p[0];
391 Nodes.IntPoint(3).x = p[1];
392 Nodes.IntPoint(3).y = p[1];
393 Nodes.IntPoint(4).x = 1. - 2. * p[1];
394 Nodes.IntPoint(4).y = p[1];
395 Nodes.IntPoint(5).x = p[1];
396 Nodes.IntPoint(5).y = 1. - 2. * p[1];
397
398 for (int i = 0; i < 6; i++)
399 {
400 const real_t x = Nodes.IntPoint(i).x, y = Nodes.IntPoint(i).y;
401 A(0,i) = 1.;
402 A(1,i) = x;
403 A(2,i) = y;
404 A(3,i) = x * x;
405 A(4,i) = x * y;
406 A(5,i) = y * y;
407 }
408
409 A.Invert();
410}
411
413 Vector &shape) const
414{
415 const real_t x = ip.x, y = ip.y;
416 pol(0) = 1.;
417 pol(1) = x;
418 pol(2) = y;
419 pol(3) = x * x;
420 pol(4) = x * y;
421 pol(5) = y * y;
422
423 A.Mult(pol, shape);
424}
425
427 DenseMatrix &dshape) const
428{
429 const real_t x = ip.x, y = ip.y;
430 D(0,0) = 0.; D(0,1) = 0.;
431 D(1,0) = 1.; D(1,1) = 0.;
432 D(2,0) = 0.; D(2,1) = 1.;
433 D(3,0) = 2. * x; D(3,1) = 0.;
434 D(4,0) = y; D(4,1) = x;
435 D(5,0) = 0.; D(5,1) = 2. * y;
436
437 Mult(A, D, dshape);
438}
439
440
442 : NodalFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
443{
444 Nodes.IntPoint(0).x = 0.0;
445 Nodes.IntPoint(0).y = 0.0;
446 Nodes.IntPoint(1).x = 1.0;
447 Nodes.IntPoint(1).y = 0.0;
448 Nodes.IntPoint(2).x = 1.0;
449 Nodes.IntPoint(2).y = 1.0;
450 Nodes.IntPoint(3).x = 0.0;
451 Nodes.IntPoint(3).y = 1.0;
452 Nodes.IntPoint(4).x = 0.5;
453 Nodes.IntPoint(4).y = 0.0;
454 Nodes.IntPoint(5).x = 1.0;
455 Nodes.IntPoint(5).y = 0.5;
456 Nodes.IntPoint(6).x = 0.5;
457 Nodes.IntPoint(6).y = 1.0;
458 Nodes.IntPoint(7).x = 0.0;
459 Nodes.IntPoint(7).y = 0.5;
460 Nodes.IntPoint(8).x = 0.5;
461 Nodes.IntPoint(8).y = 0.5;
462}
463
465 Vector &shape) const
466{
467 real_t x = ip.x, y = ip.y;
468 real_t l1x, l2x, l3x, l1y, l2y, l3y;
469
470 l1x = (x - 1.) * (2. * x - 1);
471 l2x = 4. * x * (1. - x);
472 l3x = x * (2. * x - 1.);
473 l1y = (y - 1.) * (2. * y - 1);
474 l2y = 4. * y * (1. - y);
475 l3y = y * (2. * y - 1.);
476
477 shape(0) = l1x * l1y;
478 shape(4) = l2x * l1y;
479 shape(1) = l3x * l1y;
480 shape(7) = l1x * l2y;
481 shape(8) = l2x * l2y;
482 shape(5) = l3x * l2y;
483 shape(3) = l1x * l3y;
484 shape(6) = l2x * l3y;
485 shape(2) = l3x * l3y;
486}
487
489 DenseMatrix &dshape) const
490{
491 real_t x = ip.x, y = ip.y;
492 real_t l1x, l2x, l3x, l1y, l2y, l3y;
493 real_t d1x, d2x, d3x, d1y, d2y, d3y;
494
495 l1x = (x - 1.) * (2. * x - 1);
496 l2x = 4. * x * (1. - x);
497 l3x = x * (2. * x - 1.);
498 l1y = (y - 1.) * (2. * y - 1);
499 l2y = 4. * y * (1. - y);
500 l3y = y * (2. * y - 1.);
501
502 d1x = 4. * x - 3.;
503 d2x = 4. - 8. * x;
504 d3x = 4. * x - 1.;
505 d1y = 4. * y - 3.;
506 d2y = 4. - 8. * y;
507 d3y = 4. * y - 1.;
508
509 dshape(0,0) = d1x * l1y;
510 dshape(0,1) = l1x * d1y;
511
512 dshape(4,0) = d2x * l1y;
513 dshape(4,1) = l2x * d1y;
514
515 dshape(1,0) = d3x * l1y;
516 dshape(1,1) = l3x * d1y;
517
518 dshape(7,0) = d1x * l2y;
519 dshape(7,1) = l1x * d2y;
520
521 dshape(8,0) = d2x * l2y;
522 dshape(8,1) = l2x * d2y;
523
524 dshape(5,0) = d3x * l2y;
525 dshape(5,1) = l3x * d2y;
526
527 dshape(3,0) = d1x * l3y;
528 dshape(3,1) = l1x * d3y;
529
530 dshape(6,0) = d2x * l3y;
531 dshape(6,1) = l2x * d3y;
532
533 dshape(2,0) = d3x * l3y;
534 dshape(2,1) = l3x * d3y;
535}
536
537void BiQuad2DFiniteElement::ProjectDelta(int vertex, Vector &dofs) const
538{
539#if 0
540 dofs = 1.;
541#else
542 dofs = 0.;
543 dofs(vertex) = 1.;
544 switch (vertex)
545 {
546 case 0: dofs(4) = 0.25; dofs(7) = 0.25; break;
547 case 1: dofs(4) = 0.25; dofs(5) = 0.25; break;
548 case 2: dofs(5) = 0.25; dofs(6) = 0.25; break;
549 case 3: dofs(6) = 0.25; dofs(7) = 0.25; break;
550 }
551 dofs(8) = 1./16.;
552#endif
553}
554
555
557 : NodalFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
558{
559 const real_t p1 = 0.5*(1.-sqrt(3./5.));
560
561 Nodes.IntPoint(0).x = p1;
562 Nodes.IntPoint(0).y = p1;
563 Nodes.IntPoint(4).x = 0.5;
564 Nodes.IntPoint(4).y = p1;
565 Nodes.IntPoint(1).x = 1.-p1;
566 Nodes.IntPoint(1).y = p1;
567 Nodes.IntPoint(7).x = p1;
568 Nodes.IntPoint(7).y = 0.5;
569 Nodes.IntPoint(8).x = 0.5;
570 Nodes.IntPoint(8).y = 0.5;
571 Nodes.IntPoint(5).x = 1.-p1;
572 Nodes.IntPoint(5).y = 0.5;
573 Nodes.IntPoint(3).x = p1;
574 Nodes.IntPoint(3).y = 1.-p1;
575 Nodes.IntPoint(6).x = 0.5;
576 Nodes.IntPoint(6).y = 1.-p1;
577 Nodes.IntPoint(2).x = 1.-p1;
578 Nodes.IntPoint(2).y = 1.-p1;
579}
580
582 Vector &shape) const
583{
584 const real_t a = sqrt(5./3.);
585 const real_t p1 = 0.5*(1.-sqrt(3./5.));
586
587 real_t x = a*(ip.x-p1), y = a*(ip.y-p1);
588 real_t l1x, l2x, l3x, l1y, l2y, l3y;
589
590 l1x = (x - 1.) * (2. * x - 1);
591 l2x = 4. * x * (1. - x);
592 l3x = x * (2. * x - 1.);
593 l1y = (y - 1.) * (2. * y - 1);
594 l2y = 4. * y * (1. - y);
595 l3y = y * (2. * y - 1.);
596
597 shape(0) = l1x * l1y;
598 shape(4) = l2x * l1y;
599 shape(1) = l3x * l1y;
600 shape(7) = l1x * l2y;
601 shape(8) = l2x * l2y;
602 shape(5) = l3x * l2y;
603 shape(3) = l1x * l3y;
604 shape(6) = l2x * l3y;
605 shape(2) = l3x * l3y;
606}
607
609 DenseMatrix &dshape) const
610{
611 const real_t a = sqrt(5./3.);
612 const real_t p1 = 0.5*(1.-sqrt(3./5.));
613
614 real_t x = a*(ip.x-p1), y = a*(ip.y-p1);
615 real_t l1x, l2x, l3x, l1y, l2y, l3y;
616 real_t d1x, d2x, d3x, d1y, d2y, d3y;
617
618 l1x = (x - 1.) * (2. * x - 1);
619 l2x = 4. * x * (1. - x);
620 l3x = x * (2. * x - 1.);
621 l1y = (y - 1.) * (2. * y - 1);
622 l2y = 4. * y * (1. - y);
623 l3y = y * (2. * y - 1.);
624
625 d1x = a * (4. * x - 3.);
626 d2x = a * (4. - 8. * x);
627 d3x = a * (4. * x - 1.);
628 d1y = a * (4. * y - 3.);
629 d2y = a * (4. - 8. * y);
630 d3y = a * (4. * y - 1.);
631
632 dshape(0,0) = d1x * l1y;
633 dshape(0,1) = l1x * d1y;
634
635 dshape(4,0) = d2x * l1y;
636 dshape(4,1) = l2x * d1y;
637
638 dshape(1,0) = d3x * l1y;
639 dshape(1,1) = l3x * d1y;
640
641 dshape(7,0) = d1x * l2y;
642 dshape(7,1) = l1x * d2y;
643
644 dshape(8,0) = d2x * l2y;
645 dshape(8,1) = l2x * d2y;
646
647 dshape(5,0) = d3x * l2y;
648 dshape(5,1) = l3x * d2y;
649
650 dshape(3,0) = d1x * l3y;
651 dshape(3,1) = l1x * d3y;
652
653 dshape(6,0) = d2x * l3y;
654 dshape(6,1) = l2x * d3y;
655
656 dshape(2,0) = d3x * l3y;
657 dshape(2,1) = l3x * d3y;
658}
659
661 : NodalFiniteElement (2, Geometry::SQUARE, 16, 3, FunctionSpace::Qk)
662{
663 Nodes.IntPoint(0).x = 0.;
664 Nodes.IntPoint(0).y = 0.;
665 Nodes.IntPoint(1).x = 1.;
666 Nodes.IntPoint(1).y = 0.;
667 Nodes.IntPoint(2).x = 1.;
668 Nodes.IntPoint(2).y = 1.;
669 Nodes.IntPoint(3).x = 0.;
670 Nodes.IntPoint(3).y = 1.;
671 Nodes.IntPoint(4).x = 1./3.;
672 Nodes.IntPoint(4).y = 0.;
673 Nodes.IntPoint(5).x = 2./3.;
674 Nodes.IntPoint(5).y = 0.;
675 Nodes.IntPoint(6).x = 1.;
676 Nodes.IntPoint(6).y = 1./3.;
677 Nodes.IntPoint(7).x = 1.;
678 Nodes.IntPoint(7).y = 2./3.;
679 Nodes.IntPoint(8).x = 2./3.;
680 Nodes.IntPoint(8).y = 1.;
681 Nodes.IntPoint(9).x = 1./3.;
682 Nodes.IntPoint(9).y = 1.;
683 Nodes.IntPoint(10).x = 0.;
684 Nodes.IntPoint(10).y = 2./3.;
685 Nodes.IntPoint(11).x = 0.;
686 Nodes.IntPoint(11).y = 1./3.;
687 Nodes.IntPoint(12).x = 1./3.;
688 Nodes.IntPoint(12).y = 1./3.;
689 Nodes.IntPoint(13).x = 2./3.;
690 Nodes.IntPoint(13).y = 1./3.;
691 Nodes.IntPoint(14).x = 1./3.;
692 Nodes.IntPoint(14).y = 2./3.;
693 Nodes.IntPoint(15).x = 2./3.;
694 Nodes.IntPoint(15).y = 2./3.;
695}
696
698 const IntegrationPoint &ip, Vector &shape) const
699{
700 real_t x = ip.x, y = ip.y;
701
702 real_t w1x, w2x, w3x, w1y, w2y, w3y;
703 real_t l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
704
705 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
706 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
707
708 l0x = (- 4.5) * w1x * w2x * w3x;
709 l1x = ( 13.5) * x * w2x * w3x;
710 l2x = (-13.5) * x * w1x * w3x;
711 l3x = ( 4.5) * x * w1x * w2x;
712
713 l0y = (- 4.5) * w1y * w2y * w3y;
714 l1y = ( 13.5) * y * w2y * w3y;
715 l2y = (-13.5) * y * w1y * w3y;
716 l3y = ( 4.5) * y * w1y * w2y;
717
718 shape(0) = l0x * l0y;
719 shape(1) = l3x * l0y;
720 shape(2) = l3x * l3y;
721 shape(3) = l0x * l3y;
722 shape(4) = l1x * l0y;
723 shape(5) = l2x * l0y;
724 shape(6) = l3x * l1y;
725 shape(7) = l3x * l2y;
726 shape(8) = l2x * l3y;
727 shape(9) = l1x * l3y;
728 shape(10) = l0x * l2y;
729 shape(11) = l0x * l1y;
730 shape(12) = l1x * l1y;
731 shape(13) = l2x * l1y;
732 shape(14) = l1x * l2y;
733 shape(15) = l2x * l2y;
734}
735
737 const IntegrationPoint &ip, DenseMatrix &dshape) const
738{
739 real_t x = ip.x, y = ip.y;
740
741 real_t w1x, w2x, w3x, w1y, w2y, w3y;
742 real_t l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
743 real_t d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
744
745 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
746 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
747
748 l0x = (- 4.5) * w1x * w2x * w3x;
749 l1x = ( 13.5) * x * w2x * w3x;
750 l2x = (-13.5) * x * w1x * w3x;
751 l3x = ( 4.5) * x * w1x * w2x;
752
753 l0y = (- 4.5) * w1y * w2y * w3y;
754 l1y = ( 13.5) * y * w2y * w3y;
755 l2y = (-13.5) * y * w1y * w3y;
756 l3y = ( 4.5) * y * w1y * w2y;
757
758 d0x = -5.5 + ( 18. - 13.5 * x) * x;
759 d1x = 9. + (-45. + 40.5 * x) * x;
760 d2x = -4.5 + ( 36. - 40.5 * x) * x;
761 d3x = 1. + (- 9. + 13.5 * x) * x;
762
763 d0y = -5.5 + ( 18. - 13.5 * y) * y;
764 d1y = 9. + (-45. + 40.5 * y) * y;
765 d2y = -4.5 + ( 36. - 40.5 * y) * y;
766 d3y = 1. + (- 9. + 13.5 * y) * y;
767
768 dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
769 dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
770 dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
771 dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
772 dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
773 dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
774 dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
775 dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
776 dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
777 dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
778 dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
779 dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
780 dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
781 dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
782 dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
783 dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
784}
785
787 const IntegrationPoint &ip, DenseMatrix &h) const
788{
789 real_t x = ip.x, y = ip.y;
790
791 real_t w1x, w2x, w3x, w1y, w2y, w3y;
792 real_t l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
793 real_t d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
794 real_t h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
795
796 w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
797 w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
798
799 l0x = (- 4.5) * w1x * w2x * w3x;
800 l1x = ( 13.5) * x * w2x * w3x;
801 l2x = (-13.5) * x * w1x * w3x;
802 l3x = ( 4.5) * x * w1x * w2x;
803
804 l0y = (- 4.5) * w1y * w2y * w3y;
805 l1y = ( 13.5) * y * w2y * w3y;
806 l2y = (-13.5) * y * w1y * w3y;
807 l3y = ( 4.5) * y * w1y * w2y;
808
809 d0x = -5.5 + ( 18. - 13.5 * x) * x;
810 d1x = 9. + (-45. + 40.5 * x) * x;
811 d2x = -4.5 + ( 36. - 40.5 * x) * x;
812 d3x = 1. + (- 9. + 13.5 * x) * x;
813
814 d0y = -5.5 + ( 18. - 13.5 * y) * y;
815 d1y = 9. + (-45. + 40.5 * y) * y;
816 d2y = -4.5 + ( 36. - 40.5 * y) * y;
817 d3y = 1. + (- 9. + 13.5 * y) * y;
818
819 h0x = -27. * x + 18.;
820 h1x = 81. * x - 45.;
821 h2x = -81. * x + 36.;
822 h3x = 27. * x - 9.;
823
824 h0y = -27. * y + 18.;
825 h1y = 81. * y - 45.;
826 h2y = -81. * y + 36.;
827 h3y = 27. * y - 9.;
828
829 h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
830 h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
831 h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
832 h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
833 h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
834 h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
835 h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
836 h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
837 h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
838 h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
839 h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
840 h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
841 h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
842 h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
843 h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
844 h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
845}
846
847
849 : NodalFiniteElement(1, Geometry::SEGMENT, 4, 3)
850{
851 Nodes.IntPoint(0).x = 0.0;
852 Nodes.IntPoint(1).x = 1.0;
853 Nodes.IntPoint(2).x = 0.33333333333333333333;
854 Nodes.IntPoint(3).x = 0.66666666666666666667;
855}
856
858 Vector &shape) const
859{
860 real_t x = ip.x;
861 real_t l1 = x,
862 l2 = (1.0-x),
863 l3 = (0.33333333333333333333-x),
864 l4 = (0.66666666666666666667-x);
865
866 shape(0) = 4.5 * l2 * l3 * l4;
867 shape(1) = 4.5 * l1 * l3 * l4;
868 shape(2) = 13.5 * l1 * l2 * l4;
869 shape(3) = -13.5 * l1 * l2 * l3;
870}
871
873 DenseMatrix &dshape) const
874{
875 real_t x = ip.x;
876
877 dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
878 dshape(1,0) = 1. - x * (9. - 13.5 * x);
879 dshape(2,0) = 9. - x * (45. - 40.5 * x);
880 dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
881}
882
883
885 : NodalFiniteElement(2, Geometry::TRIANGLE, 10, 3)
886{
887 Nodes.IntPoint(0).x = 0.0;
888 Nodes.IntPoint(0).y = 0.0;
889 Nodes.IntPoint(1).x = 1.0;
890 Nodes.IntPoint(1).y = 0.0;
891 Nodes.IntPoint(2).x = 0.0;
892 Nodes.IntPoint(2).y = 1.0;
893 Nodes.IntPoint(3).x = 0.33333333333333333333;
894 Nodes.IntPoint(3).y = 0.0;
895 Nodes.IntPoint(4).x = 0.66666666666666666667;
896 Nodes.IntPoint(4).y = 0.0;
897 Nodes.IntPoint(5).x = 0.66666666666666666667;
898 Nodes.IntPoint(5).y = 0.33333333333333333333;
899 Nodes.IntPoint(6).x = 0.33333333333333333333;
900 Nodes.IntPoint(6).y = 0.66666666666666666667;
901 Nodes.IntPoint(7).x = 0.0;
902 Nodes.IntPoint(7).y = 0.66666666666666666667;
903 Nodes.IntPoint(8).x = 0.0;
904 Nodes.IntPoint(8).y = 0.33333333333333333333;
905 Nodes.IntPoint(9).x = 0.33333333333333333333;
906 Nodes.IntPoint(9).y = 0.33333333333333333333;
907}
908
910 Vector &shape) const
911{
912 real_t x = ip.x, y = ip.y;
913 real_t l1 = (-1. + x + y),
914 lx = (-1. + 3.*x),
915 ly = (-1. + 3.*y);
916
917 shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
918 shape(1) = 0.5*x*(lx - 1.)*lx;
919 shape(2) = 0.5*y*(-1. + ly)*ly;
920 shape(3) = 4.5*x*l1*(3.*l1 + 1.);
921 shape(4) = -4.5*x*lx*l1;
922 shape(5) = 4.5*x*lx*y;
923 shape(6) = 4.5*x*y*ly;
924 shape(7) = -4.5*y*l1*ly;
925 shape(8) = 4.5*y*l1*(1. + 3.*l1);
926 shape(9) = -27.*x*y*l1;
927}
928
930 DenseMatrix &dshape) const
931{
932 real_t x = ip.x, y = ip.y;
933
934 dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
935 dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
936 dshape(2,0) = 0.;
937 dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
938 dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
939 dshape(5,0) = 4.5*(-1. + 6.*x)*y;
940 dshape(6,0) = 4.5*y*(-1. + 3.*y);
941 dshape(7,0) = 4.5*(1. - 3.*y)*y;
942 dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
943 dshape(9,0) = -27.*y*(-1. + 2.*x + y);
944
945 dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
946 dshape(1,1) = 0.;
947 dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
948 dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
949 dshape(4,1) = 4.5*(1. - 3.*x)*x;
950 dshape(5,1) = 4.5*x*(-1. + 3.*x);
951 dshape(6,1) = 4.5*x*(-1. + 6.*y);
952 dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
953 dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
954 dshape(9,1) = -27.*x*(-1. + x + 2.*y);
955}
956
958 DenseMatrix &h) const
959{
960 real_t x = ip.x, y = ip.y;
961
962 h(0,0) = 18.-27.*(x+y);
963 h(0,1) = 18.-27.*(x+y);
964 h(0,2) = 18.-27.*(x+y);
965
966 h(1,0) = -9.+27.*x;
967 h(1,1) = 0.;
968 h(1,2) = 0.;
969
970 h(2,0) = 0.;
971 h(2,1) = 0.;
972 h(2,2) = -9.+27.*y;
973
974 h(3,0) = -45.+81.*x+54.*y;
975 h(3,1) = -22.5+54.*x+27.*y;
976 h(3,2) = 27.*x;
977
978 h(4,0) = 36.-81.*x-27.*y;
979 h(4,1) = 4.5-27.*x;
980 h(4,2) = 0.;
981
982 h(5,0) = 27.*y;
983 h(5,1) = -4.5+27.*x;
984 h(5,2) = 0.;
985
986 h(6,0) = 0.;
987 h(6,1) = -4.5+27.*y;
988 h(6,2) = 27.*x;
989
990 h(7,0) = 0.;
991 h(7,1) = 4.5-27.*y;
992 h(7,2) = 36.-27.*x-81.*y;
993
994 h(8,0) = 27.*y;
995 h(8,1) = -22.5+27.*x+54.*y;
996 h(8,2) = -45.+54.*x+81.*y;
997
998 h(9,0) = -54.*y;
999 h(9,1) = 27.-54.*(x+y);
1000 h(9,2) = -54.*x;
1001}
1002
1003
1005 : NodalFiniteElement(3, Geometry::TETRAHEDRON, 20, 3)
1006{
1007 Nodes.IntPoint(0).x = 0;
1008 Nodes.IntPoint(0).y = 0;
1009 Nodes.IntPoint(0).z = 0;
1010 Nodes.IntPoint(1).x = 1.;
1011 Nodes.IntPoint(1).y = 0;
1012 Nodes.IntPoint(1).z = 0;
1013 Nodes.IntPoint(2).x = 0;
1014 Nodes.IntPoint(2).y = 1.;
1015 Nodes.IntPoint(2).z = 0;
1016 Nodes.IntPoint(3).x = 0;
1017 Nodes.IntPoint(3).y = 0;
1018 Nodes.IntPoint(3).z = 1.;
1019 Nodes.IntPoint(4).x = 0.3333333333333333333333333333;
1020 Nodes.IntPoint(4).y = 0;
1021 Nodes.IntPoint(4).z = 0;
1022 Nodes.IntPoint(5).x = 0.6666666666666666666666666667;
1023 Nodes.IntPoint(5).y = 0;
1024 Nodes.IntPoint(5).z = 0;
1025 Nodes.IntPoint(6).x = 0;
1026 Nodes.IntPoint(6).y = 0.3333333333333333333333333333;
1027 Nodes.IntPoint(6).z = 0;
1028 Nodes.IntPoint(7).x = 0;
1029 Nodes.IntPoint(7).y = 0.6666666666666666666666666667;
1030 Nodes.IntPoint(7).z = 0;
1031 Nodes.IntPoint(8).x = 0;
1032 Nodes.IntPoint(8).y = 0;
1033 Nodes.IntPoint(8).z = 0.3333333333333333333333333333;
1034 Nodes.IntPoint(9).x = 0;
1035 Nodes.IntPoint(9).y = 0;
1036 Nodes.IntPoint(9).z = 0.6666666666666666666666666667;
1037 Nodes.IntPoint(10).x = 0.6666666666666666666666666667;
1038 Nodes.IntPoint(10).y = 0.3333333333333333333333333333;
1039 Nodes.IntPoint(10).z = 0;
1040 Nodes.IntPoint(11).x = 0.3333333333333333333333333333;
1041 Nodes.IntPoint(11).y = 0.6666666666666666666666666667;
1042 Nodes.IntPoint(11).z = 0;
1043 Nodes.IntPoint(12).x = 0.6666666666666666666666666667;
1044 Nodes.IntPoint(12).y = 0;
1045 Nodes.IntPoint(12).z = 0.3333333333333333333333333333;
1046 Nodes.IntPoint(13).x = 0.3333333333333333333333333333;
1047 Nodes.IntPoint(13).y = 0;
1048 Nodes.IntPoint(13).z = 0.6666666666666666666666666667;
1049 Nodes.IntPoint(14).x = 0;
1050 Nodes.IntPoint(14).y = 0.6666666666666666666666666667;
1051 Nodes.IntPoint(14).z = 0.3333333333333333333333333333;
1052 Nodes.IntPoint(15).x = 0;
1053 Nodes.IntPoint(15).y = 0.3333333333333333333333333333;
1054 Nodes.IntPoint(15).z = 0.6666666666666666666666666667;
1055 Nodes.IntPoint(16).x = 0.3333333333333333333333333333;
1056 Nodes.IntPoint(16).y = 0.3333333333333333333333333333;
1057 Nodes.IntPoint(16).z = 0.3333333333333333333333333333;
1058 Nodes.IntPoint(17).x = 0;
1059 Nodes.IntPoint(17).y = 0.3333333333333333333333333333;
1060 Nodes.IntPoint(17).z = 0.3333333333333333333333333333;
1061 Nodes.IntPoint(18).x = 0.3333333333333333333333333333;
1062 Nodes.IntPoint(18).y = 0;
1063 Nodes.IntPoint(18).z = 0.3333333333333333333333333333;
1064 Nodes.IntPoint(19).x = 0.3333333333333333333333333333;
1065 Nodes.IntPoint(19).y = 0.3333333333333333333333333333;
1066 Nodes.IntPoint(19).z = 0;
1067}
1068
1070 Vector &shape) const
1071{
1072 real_t x = ip.x, y = ip.y, z = ip.z;
1073
1074 shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
1075 (-1 + 3*x + 3*y + 3*z))/2.;
1076 shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1077 shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
1078 shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
1079 shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1080 shape(19) = -27*x*y*(-1 + x + y + z);
1081 shape(10) = (9*x*(-1 + 3*x)*y)/2.;
1082 shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
1083 shape(11) = (9*x*y*(-1 + 3*y))/2.;
1084 shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
1085 shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
1086 shape(18) = -27*x*z*(-1 + x + y + z);
1087 shape(12) = (9*x*(-1 + 3*x)*z)/2.;
1088 shape(17) = -27*y*z*(-1 + x + y + z);
1089 shape(16) = 27*x*y*z;
1090 shape(14) = (9*y*(-1 + 3*y)*z)/2.;
1091 shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
1092 shape(13) = (9*x*z*(-1 + 3*z))/2.;
1093 shape(15) = (9*y*z*(-1 + 3*z))/2.;
1094 shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
1095}
1096
1098 DenseMatrix &dshape) const
1099{
1100 real_t x = ip.x, y = ip.y, z = ip.z;
1101
1102 dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1103 x*(-4 + 6*y + 6*z)))/2.;
1104 dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1105 x*(-4 + 6*y + 6*z)))/2.;
1106 dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
1107 x*(-4 + 6*y + 6*z)))/2.;
1108 dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
1109 2*x*(-5 + 6*y + 6*z)))/2.;
1110 dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
1111 dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
1112 dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
1113 dshape(5,1) = (9*(1 - 3*x)*x)/2.;
1114 dshape(5,2) = (9*(1 - 3*x)*x)/2.;
1115 dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
1116 dshape(1,1) = 0;
1117 dshape(1,2) = 0;
1118 dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
1119 dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
1120 x*(-5 + 12*y + 6*z)))/2.;
1121 dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
1122 dshape(19,0) = -27*y*(-1 + 2*x + y + z);
1123 dshape(19,1) = -27*x*(-1 + x + 2*y + z);
1124 dshape(19,2) = -27*x*y;
1125 dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
1126 dshape(10,1) = (9*x*(-1 + 3*x))/2.;
1127 dshape(10,2) = 0;
1128 dshape(7,0) = (9*(1 - 3*y)*y)/2.;
1129 dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
1130 dshape(7,2) = (9*(1 - 3*y)*y)/2.;
1131 dshape(11,0) = (9*y*(-1 + 3*y))/2.;
1132 dshape(11,1) = (9*x*(-1 + 6*y))/2.;
1133 dshape(11,2) = 0;
1134 dshape(2,0) = 0;
1135 dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
1136 dshape(2,2) = 0;
1137 dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
1138 dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
1139 dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
1140 x*(-5 + 6*y + 12*z)))/2.;
1141 dshape(18,0) = -27*z*(-1 + 2*x + y + z);
1142 dshape(18,1) = -27*x*z;
1143 dshape(18,2) = -27*x*(-1 + x + y + 2*z);
1144 dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
1145 dshape(12,1) = 0;
1146 dshape(12,2) = (9*x*(-1 + 3*x))/2.;
1147 dshape(17,0) = -27*y*z;
1148 dshape(17,1) = -27*z*(-1 + x + 2*y + z);
1149 dshape(17,2) = -27*y*(-1 + x + y + 2*z);
1150 dshape(16,0) = 27*y*z;
1151 dshape(16,1) = 27*x*z;
1152 dshape(16,2) = 27*x*y;
1153 dshape(14,0) = 0;
1154 dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
1155 dshape(14,2) = (9*y*(-1 + 3*y))/2.;
1156 dshape(9,0) = (9*(1 - 3*z)*z)/2.;
1157 dshape(9,1) = (9*(1 - 3*z)*z)/2.;
1158 dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
1159 dshape(13,0) = (9*z*(-1 + 3*z))/2.;
1160 dshape(13,1) = 0;
1161 dshape(13,2) = (9*x*(-1 + 6*z))/2.;
1162 dshape(15,0) = 0;
1163 dshape(15,1) = (9*z*(-1 + 3*z))/2.;
1164 dshape(15,2) = (9*y*(-1 + 6*z))/2.;
1165 dshape(3,0) = 0;
1166 dshape(3,1) = 0;
1167 dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
1168}
1169
1170
1172 : NodalFiniteElement(2, Geometry::TRIANGLE, 1, 0)
1173{
1174 Nodes.IntPoint(0).x = 0.333333333333333333;
1175 Nodes.IntPoint(0).y = 0.333333333333333333;
1176}
1177
1179 Vector &shape) const
1180{
1181 shape(0) = 1.0;
1182}
1183
1185 DenseMatrix &dshape) const
1186{
1187 dshape(0,0) = 0.0;
1188 dshape(0,1) = 0.0;
1189}
1190
1191
1193 : NodalFiniteElement(2, Geometry::SQUARE, 1, 0, FunctionSpace::Qk)
1194{
1195 Nodes.IntPoint(0).x = 0.5;
1196 Nodes.IntPoint(0).y = 0.5;
1197}
1198
1200 Vector &shape) const
1201{
1202 shape(0) = 1.0;
1203}
1204
1206 DenseMatrix &dshape) const
1207{
1208 dshape(0,0) = 0.0;
1209 dshape(0,1) = 0.0;
1210}
1211
1212
1214 : NodalFiniteElement(3, Geometry::TETRAHEDRON, 4, 1)
1215{
1216 Nodes.IntPoint(0).x = 0.0;
1217 Nodes.IntPoint(0).y = 0.0;
1218 Nodes.IntPoint(0).z = 0.0;
1219 Nodes.IntPoint(1).x = 1.0;
1220 Nodes.IntPoint(1).y = 0.0;
1221 Nodes.IntPoint(1).z = 0.0;
1222 Nodes.IntPoint(2).x = 0.0;
1223 Nodes.IntPoint(2).y = 1.0;
1224 Nodes.IntPoint(2).z = 0.0;
1225 Nodes.IntPoint(3).x = 0.0;
1226 Nodes.IntPoint(3).y = 0.0;
1227 Nodes.IntPoint(3).z = 1.0;
1228}
1229
1231 Vector &shape) const
1232{
1233 shape(0) = 1. - ip.x - ip.y - ip.z;
1234 shape(1) = ip.x;
1235 shape(2) = ip.y;
1236 shape(3) = ip.z;
1237}
1238
1240 DenseMatrix &dshape) const
1241{
1242 if (dshape.Height() == 4)
1243 {
1244 real_t *A = &dshape(0,0);
1245 A[0] = -1.; A[4] = -1.; A[8] = -1.;
1246 A[1] = 1.; A[5] = 0.; A[9] = 0.;
1247 A[2] = 0.; A[6] = 1.; A[10] = 0.;
1248 A[3] = 0.; A[7] = 0.; A[11] = 1.;
1249 }
1250 else
1251 {
1252 dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
1253 dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
1254 dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
1255 dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
1256 }
1257}
1258
1259void Linear3DFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
1260const
1261{
1262 static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
1263
1264 *ndofs = 3;
1265 *dofs = face_dofs[face];
1266}
1267
1268
1269// TODO: use a FunctionSpace specific to wedges instead of Qk.
1271 : NodalFiniteElement(3, Geometry::PRISM, 6, 1, FunctionSpace::Qk)
1272{
1273 Nodes.IntPoint(0).x = 0.0;
1274 Nodes.IntPoint(0).y = 0.0;
1275 Nodes.IntPoint(0).z = 0.0;
1276 Nodes.IntPoint(1).x = 1.0;
1277 Nodes.IntPoint(1).y = 0.0;
1278 Nodes.IntPoint(1).z = 0.0;
1279 Nodes.IntPoint(2).x = 0.0;
1280 Nodes.IntPoint(2).y = 1.0;
1281 Nodes.IntPoint(2).z = 0.0;
1282 Nodes.IntPoint(3).x = 0.0;
1283 Nodes.IntPoint(3).y = 0.0;
1284 Nodes.IntPoint(3).z = 1.0;
1285 Nodes.IntPoint(4).x = 1.0;
1286 Nodes.IntPoint(4).y = 0.0;
1287 Nodes.IntPoint(4).z = 1.0;
1288 Nodes.IntPoint(5).x = 0.0;
1289 Nodes.IntPoint(5).y = 1.0;
1290 Nodes.IntPoint(5).z = 1.0;
1291}
1292
1294 Vector &shape) const
1295{
1296 shape(0) = (1. - ip.x - ip.y) * (1. - ip.z);
1297 shape(1) = ip.x * (1. - ip.z);
1298 shape(2) = ip.y * (1. - ip.z);
1299 shape(3) = (1. - ip.x - ip.y) * ip.z;
1300 shape(4) = ip.x * ip.z;
1301 shape(5) = ip.y * ip.z;
1302}
1303
1305 DenseMatrix &dshape) const
1306{
1307 dshape(0,0) = -1. + ip.z;
1308 dshape(0,1) = -1. + ip.z;
1309 dshape(0,2) = -1. + ip.x + ip.y;
1310
1311 dshape(1,0) = 1. - ip.z;
1312 dshape(1,1) = 0.;
1313 dshape(1,2) = -ip.x;
1314
1315 dshape(2,0) = 0.;
1316 dshape(2,1) = 1. - ip.z;
1317 dshape(2,2) = -ip.y;
1318
1319 dshape(3,0) = -ip.z;
1320 dshape(3,1) = -ip.z;
1321 dshape(3,2) = 1. - ip.x - ip.y;
1322
1323 dshape(4,0) = ip.z;
1324 dshape(4,1) = 0.;
1325 dshape(4,2) = ip.x;
1326
1327 dshape(5,0) = 0.;
1328 dshape(5,1) = ip.z;
1329 dshape(5,2) = ip.y;
1330}
1331
1332void LinearWedgeFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
1333const
1334{
1335 static int face_dofs[5][4] =
1336 {{0, 2, 1, -1}, {3, 4, 5, -1}, {0, 1, 4, 3}, {1, 2, 5, 4}, {2, 0, 3, 5}};
1337
1338 *ndofs = (face < 2) ? 3 : 4;
1339 *dofs = face_dofs[face];
1340}
1341
1342
1344 : NodalFiniteElement(3, Geometry::PYRAMID, 5, 1)
1345{
1346 Nodes.IntPoint(0).x = 0.0;
1347 Nodes.IntPoint(0).y = 0.0;
1348 Nodes.IntPoint(0).z = 0.0;
1349 Nodes.IntPoint(1).x = 1.0;
1350 Nodes.IntPoint(1).y = 0.0;
1351 Nodes.IntPoint(1).z = 0.0;
1352 Nodes.IntPoint(2).x = 1.0;
1353 Nodes.IntPoint(2).y = 1.0;
1354 Nodes.IntPoint(2).z = 0.0;
1355 Nodes.IntPoint(3).x = 0.0;
1356 Nodes.IntPoint(3).y = 1.0;
1357 Nodes.IntPoint(3).z = 0.0;
1358 Nodes.IntPoint(4).x = 0.0;
1359 Nodes.IntPoint(4).y = 0.0;
1360 Nodes.IntPoint(4).z = 1.0;
1361}
1362
1364 Vector &shape) const
1365{
1366 real_t x = ip.x, y = ip.y, z = ip.z;
1367 real_t ox = 1.-x-z, oy = 1.-y-z, oz = 1.-z;
1368
1369 real_t tol = 1e-6;
1370
1371 if (oz <= tol)
1372 {
1373 // We must return the limit of the basis functions as z->1. In order to
1374 // remain inside the pyramid in this limit the x and y coordinates must
1375 // be approaching 0. The resulting limiting basis function values are:
1376 shape(0) = 0.;
1377 shape(1) = 0.;
1378 shape(2) = 0.;
1379 shape(3) = 0.;
1380 shape(4) = 1.;
1381 return;
1382 }
1383
1384 real_t ozi = 1. / oz;
1385
1386 shape(0) = ox * oy * ozi;
1387 shape(1) = x * oy * ozi;
1388 shape(2) = x * y * ozi;
1389 shape(3) = ox * y * ozi;
1390 shape(4) = z;
1391}
1392
1394 DenseMatrix &dshape) const
1395{
1396 real_t x = ip.x, y = ip.y, z = ip.z;
1397 real_t ox = 1.-x-z, oy = 1.-y-z, oz = 1.-z;
1398
1399 real_t tol = 1e-6;
1400
1401 if (oz <= tol)
1402 {
1403 // At the apex of the pyramid the gradients of the basis functions are
1404 // multivalued and depend on the direction from which the limit is taken.
1405 // The following values correspond to the average of the gradients taken
1406 // over all possible directions approaching the apex of the pyramid from
1407 // within its interior.
1408 dshape(0,0) = - 0.5;
1409 dshape(0,1) = - 0.5;
1410 dshape(0,2) = - 0.75;
1411
1412 dshape(1,0) = 0.5;
1413 dshape(1,1) = - 0.5;
1414 dshape(1,2) = - 0.25;
1415
1416 dshape(2,0) = 0.5;
1417 dshape(2,1) = 0.5;
1418 dshape(2,2) = 0.25;
1419
1420 dshape(3,0) = - 0.5;
1421 dshape(3,1) = 0.5;
1422 dshape(3,2) = - 0.25;
1423
1424 dshape(4,0) = 0.;
1425 dshape(4,1) = 0.;
1426 dshape(4,2) = 1.;
1427
1428 return;
1429 }
1430
1431 real_t ozi = 1. / oz;
1432
1433 dshape(0,0) = - oy * ozi;
1434 dshape(0,1) = - ox * ozi;
1435 dshape(0,2) = x * y * ozi * ozi - 1.;
1436
1437 dshape(1,0) = oy * ozi;
1438 dshape(1,1) = - x * ozi;
1439 dshape(1,2) = - x * y * ozi * ozi;
1440
1441 dshape(2,0) = y * ozi;
1442 dshape(2,1) = x * ozi;
1443 dshape(2,2) = x * y * ozi * ozi;
1444
1445 dshape(3,0) = - y * ozi;
1446 dshape(3,1) = ox * ozi;
1447 dshape(3,2) = - x * y * ozi * ozi;
1448
1449 dshape(4,0) = 0.;
1450 dshape(4,1) = 0.;
1451 dshape(4,2) = 1.;
1452}
1453
1454void LinearPyramidFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
1455const
1456{
1457 static int face_dofs[5][4] =
1458 {{3, 2, 1, 0}, {0, 1, 4, -1}, {1, 2, 4, -1}, {2, 3, 4, -1}, {3, 0, 4, -1}};
1459
1460 *ndofs = (face < 1) ? 4 : 3;
1461 *dofs = face_dofs[face];
1462}
1463
1464
1466 : NodalFiniteElement(3, Geometry::TETRAHEDRON, 10, 2)
1467{
1468 Nodes.IntPoint(0).x = 0.0;
1469 Nodes.IntPoint(0).y = 0.0;
1470 Nodes.IntPoint(0).z = 0.0;
1471 Nodes.IntPoint(1).x = 1.0;
1472 Nodes.IntPoint(1).y = 0.0;
1473 Nodes.IntPoint(1).z = 0.0;
1474 Nodes.IntPoint(2).x = 0.0;
1475 Nodes.IntPoint(2).y = 1.0;
1476 Nodes.IntPoint(2).z = 0.0;
1477 Nodes.IntPoint(3).x = 0.0;
1478 Nodes.IntPoint(3).y = 0.0;
1479 Nodes.IntPoint(3).z = 1.0;
1480 Nodes.IntPoint(4).x = 0.5;
1481 Nodes.IntPoint(4).y = 0.0;
1482 Nodes.IntPoint(4).z = 0.0;
1483 Nodes.IntPoint(5).x = 0.0;
1484 Nodes.IntPoint(5).y = 0.5;
1485 Nodes.IntPoint(5).z = 0.0;
1486 Nodes.IntPoint(6).x = 0.0;
1487 Nodes.IntPoint(6).y = 0.0;
1488 Nodes.IntPoint(6).z = 0.5;
1489 Nodes.IntPoint(7).x = 0.5;
1490 Nodes.IntPoint(7).y = 0.5;
1491 Nodes.IntPoint(7).z = 0.0;
1492 Nodes.IntPoint(8).x = 0.5;
1493 Nodes.IntPoint(8).y = 0.0;
1494 Nodes.IntPoint(8).z = 0.5;
1495 Nodes.IntPoint(9).x = 0.0;
1496 Nodes.IntPoint(9).y = 0.5;
1497 Nodes.IntPoint(9).z = 0.5;
1498}
1499
1501 Vector &shape) const
1502{
1503 real_t L0, L1, L2, L3;
1504
1505 L0 = 1. - ip.x - ip.y - ip.z;
1506 L1 = ip.x;
1507 L2 = ip.y;
1508 L3 = ip.z;
1509
1510 shape(0) = L0 * ( 2.0 * L0 - 1.0 );
1511 shape(1) = L1 * ( 2.0 * L1 - 1.0 );
1512 shape(2) = L2 * ( 2.0 * L2 - 1.0 );
1513 shape(3) = L3 * ( 2.0 * L3 - 1.0 );
1514 shape(4) = 4.0 * L0 * L1;
1515 shape(5) = 4.0 * L0 * L2;
1516 shape(6) = 4.0 * L0 * L3;
1517 shape(7) = 4.0 * L1 * L2;
1518 shape(8) = 4.0 * L1 * L3;
1519 shape(9) = 4.0 * L2 * L3;
1520}
1521
1523 DenseMatrix &dshape) const
1524{
1525 real_t x, y, z, L0;
1526
1527 x = ip.x;
1528 y = ip.y;
1529 z = ip.z;
1530 L0 = 1.0 - x - y - z;
1531
1532 dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
1533 dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
1534 dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
1535 dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
1536 dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
1537 dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
1538 dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
1539 dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
1540 dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
1541 dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
1542}
1543
1545 : NodalFiniteElement(3, Geometry::CUBE, 8, 1, FunctionSpace::Qk)
1546{
1547 Nodes.IntPoint(0).x = 0.0;
1548 Nodes.IntPoint(0).y = 0.0;
1549 Nodes.IntPoint(0).z = 0.0;
1550
1551 Nodes.IntPoint(1).x = 1.0;
1552 Nodes.IntPoint(1).y = 0.0;
1553 Nodes.IntPoint(1).z = 0.0;
1554
1555 Nodes.IntPoint(2).x = 1.0;
1556 Nodes.IntPoint(2).y = 1.0;
1557 Nodes.IntPoint(2).z = 0.0;
1558
1559 Nodes.IntPoint(3).x = 0.0;
1560 Nodes.IntPoint(3).y = 1.0;
1561 Nodes.IntPoint(3).z = 0.0;
1562
1563 Nodes.IntPoint(4).x = 0.0;
1564 Nodes.IntPoint(4).y = 0.0;
1565 Nodes.IntPoint(4).z = 1.0;
1566
1567 Nodes.IntPoint(5).x = 1.0;
1568 Nodes.IntPoint(5).y = 0.0;
1569 Nodes.IntPoint(5).z = 1.0;
1570
1571 Nodes.IntPoint(6).x = 1.0;
1572 Nodes.IntPoint(6).y = 1.0;
1573 Nodes.IntPoint(6).z = 1.0;
1574
1575 Nodes.IntPoint(7).x = 0.0;
1576 Nodes.IntPoint(7).y = 1.0;
1577 Nodes.IntPoint(7).z = 1.0;
1578}
1579
1581 Vector &shape) const
1582{
1583 real_t x = ip.x, y = ip.y, z = ip.z;
1584 real_t ox = 1.-x, oy = 1.-y, oz = 1.-z;
1585
1586 shape(0) = ox * oy * oz;
1587 shape(1) = x * oy * oz;
1588 shape(2) = x * y * oz;
1589 shape(3) = ox * y * oz;
1590 shape(4) = ox * oy * z;
1591 shape(5) = x * oy * z;
1592 shape(6) = x * y * z;
1593 shape(7) = ox * y * z;
1594}
1595
1597 DenseMatrix &dshape) const
1598{
1599 real_t x = ip.x, y = ip.y, z = ip.z;
1600 real_t ox = 1.-x, oy = 1.-y, oz = 1.-z;
1601
1602 dshape(0,0) = - oy * oz;
1603 dshape(0,1) = - ox * oz;
1604 dshape(0,2) = - ox * oy;
1605
1606 dshape(1,0) = oy * oz;
1607 dshape(1,1) = - x * oz;
1608 dshape(1,2) = - x * oy;
1609
1610 dshape(2,0) = y * oz;
1611 dshape(2,1) = x * oz;
1612 dshape(2,2) = - x * y;
1613
1614 dshape(3,0) = - y * oz;
1615 dshape(3,1) = ox * oz;
1616 dshape(3,2) = - ox * y;
1617
1618 dshape(4,0) = - oy * z;
1619 dshape(4,1) = - ox * z;
1620 dshape(4,2) = ox * oy;
1621
1622 dshape(5,0) = oy * z;
1623 dshape(5,1) = - x * z;
1624 dshape(5,2) = x * oy;
1625
1626 dshape(6,0) = y * z;
1627 dshape(6,1) = x * z;
1628 dshape(6,2) = x * y;
1629
1630 dshape(7,0) = - y * z;
1631 dshape(7,1) = ox * z;
1632 dshape(7,2) = ox * y;
1633}
1634
1635
1637 : NodalFiniteElement(1, Geometry::SEGMENT, 1, Ord) // default Ord = 0
1638{
1639 Nodes.IntPoint(0).x = 0.5;
1640}
1641
1643 Vector &shape) const
1644{
1645 shape(0) = 1.0;
1646}
1647
1649 DenseMatrix &dshape) const
1650{
1651 dshape(0,0) = 0.0;
1652}
1653
1655 : NodalFiniteElement(2, Geometry::TRIANGLE, 3, 1)
1656{
1657 Nodes.IntPoint(0).x = 0.5;
1658 Nodes.IntPoint(0).y = 0.0;
1659 Nodes.IntPoint(1).x = 0.5;
1660 Nodes.IntPoint(1).y = 0.5;
1661 Nodes.IntPoint(2).x = 0.0;
1662 Nodes.IntPoint(2).y = 0.5;
1663}
1664
1666 Vector &shape) const
1667{
1668 shape(0) = 1.0 - 2.0 * ip.y;
1669 shape(1) = -1.0 + 2.0 * ( ip.x + ip.y );
1670 shape(2) = 1.0 - 2.0 * ip.x;
1671}
1672
1674 DenseMatrix &dshape) const
1675{
1676 dshape(0,0) = 0.0; dshape(0,1) = -2.0;
1677 dshape(1,0) = 2.0; dshape(1,1) = 2.0;
1678 dshape(2,0) = -2.0; dshape(2,1) = 0.0;
1679}
1680
1682// the FunctionSpace should be rotated (45 degrees) Q_1
1683// i.e. the span of { 1, x, y, x^2 - y^2 }
1684 : NodalFiniteElement(2, Geometry::SQUARE, 4, 2, FunctionSpace::Qk)
1685{
1686 Nodes.IntPoint(0).x = 0.5;
1687 Nodes.IntPoint(0).y = 0.0;
1688 Nodes.IntPoint(1).x = 1.0;
1689 Nodes.IntPoint(1).y = 0.5;
1690 Nodes.IntPoint(2).x = 0.5;
1691 Nodes.IntPoint(2).y = 1.0;
1692 Nodes.IntPoint(3).x = 0.0;
1693 Nodes.IntPoint(3).y = 0.5;
1694}
1695
1697 Vector &shape) const
1698{
1699 const real_t l1 = ip.x+ip.y-0.5, l2 = 1.-l1, l3 = ip.x-ip.y+0.5, l4 = 1.-l3;
1700
1701 shape(0) = l2 * l3;
1702 shape(1) = l1 * l3;
1703 shape(2) = l1 * l4;
1704 shape(3) = l2 * l4;
1705}
1706
1708 DenseMatrix &dshape) const
1709{
1710 const real_t x2 = 2.*ip.x, y2 = 2.*ip.y;
1711
1712 dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
1713 dshape(1,0) = x2; dshape(1,1) = 1. - y2;
1714 dshape(2,0) = 1. - x2; dshape(2,1) = y2;
1715 dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
1716}
1717
1718
1720 : VectorFiniteElement(2, Geometry::TRIANGLE, 3, 1, H_DIV)
1721{
1722 Nodes.IntPoint(0).x = 0.5;
1723 Nodes.IntPoint(0).y = 0.0;
1724 Nodes.IntPoint(1).x = 0.5;
1725 Nodes.IntPoint(1).y = 0.5;
1726 Nodes.IntPoint(2).x = 0.0;
1727 Nodes.IntPoint(2).y = 0.5;
1728}
1729
1731 DenseMatrix &shape) const
1732{
1733 real_t x = ip.x, y = ip.y;
1734
1735 shape(0,0) = x;
1736 shape(0,1) = y - 1.;
1737 shape(1,0) = x;
1738 shape(1,1) = y;
1739 shape(2,0) = x - 1.;
1740 shape(2,1) = y;
1741}
1742
1744 Vector &divshape) const
1745{
1746 divshape(0) = 2.;
1747 divshape(1) = 2.;
1748 divshape(2) = 2.;
1749}
1750
1751const real_t RT0TriangleFiniteElement::nk[3][2] =
1752{ {0, -1}, {1, 1}, {-1, 0} };
1753
1755 ElementTransformation &Trans, DenseMatrix &I) const
1756{
1757 int k, j;
1760#endif
1761
1762#ifdef MFEM_DEBUG
1763 for (k = 0; k < 3; k++)
1764 {
1766 for (j = 0; j < 3; j++)
1767 {
1768 real_t d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
1769 if (j == k) { d -= 1.0; }
1770 if (fabs(d) > 1.0e-12)
1771 {
1772 mfem::err << "RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
1773 " k = " << k << ", j = " << j << ", d = " << d << endl;
1774 mfem_error();
1775 }
1776 }
1777 }
1778#endif
1779
1781 ip.x = ip.y = 0.0;
1782 Trans.SetIntPoint (&ip);
1783 // Trans must be linear
1784 // set Jinv = |J| J^{-t} = adj(J)^t
1785 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
1786
1787 real_t vk[2];
1788 Vector xk (vk, 2);
1789
1790 for (k = 0; k < 3; k++)
1791 {
1792 Trans.Transform (Nodes.IntPoint (k), xk);
1793 ip.x = vk[0]; ip.y = vk[1];
1794 CalcVShape (ip, vshape);
1795 // vk = |J| J^{-t} nk
1796 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
1797 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
1798 for (j = 0; j < 3; j++)
1799 if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
1800 {
1801 I(k,j) = 0.0;
1802 }
1803 }
1804}
1805
1808 Vector &dofs) const
1809{
1810 real_t vk[2];
1811 Vector xk (vk, 2);
1812
1813 for (int k = 0; k < 3; k++)
1814 {
1815 Trans.SetIntPoint (&Nodes.IntPoint (k));
1816 // set Jinv = |J| J^{-t} = adj(J)^t
1817 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
1818
1819 vc.Eval (xk, Trans, Nodes.IntPoint (k));
1820 // xk^t |J| J^{-t} nk
1821 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
1822 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
1823 }
1824}
1825
1827 : VectorFiniteElement(2, Geometry::SQUARE, 4, 1, H_DIV,
1828 FunctionSpace::Qk)
1829{
1830 Nodes.IntPoint(0).x = 0.5;
1831 Nodes.IntPoint(0).y = 0.0;
1832 Nodes.IntPoint(1).x = 1.0;
1833 Nodes.IntPoint(1).y = 0.5;
1834 Nodes.IntPoint(2).x = 0.5;
1835 Nodes.IntPoint(2).y = 1.0;
1836 Nodes.IntPoint(3).x = 0.0;
1837 Nodes.IntPoint(3).y = 0.5;
1838}
1839
1841 DenseMatrix &shape) const
1842{
1843 real_t x = ip.x, y = ip.y;
1844
1845 shape(0,0) = 0;
1846 shape(0,1) = y - 1.;
1847 shape(1,0) = x;
1848 shape(1,1) = 0;
1849 shape(2,0) = 0;
1850 shape(2,1) = y;
1851 shape(3,0) = x - 1.;
1852 shape(3,1) = 0;
1853}
1854
1856 Vector &divshape) const
1857{
1858 divshape(0) = 1.;
1859 divshape(1) = 1.;
1860 divshape(2) = 1.;
1861 divshape(3) = 1.;
1862}
1863
1865{ {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
1866
1868 ElementTransformation &Trans, DenseMatrix &I) const
1869{
1870 int k, j;
1873#endif
1874
1875#ifdef MFEM_DEBUG
1876 for (k = 0; k < 4; k++)
1877 {
1879 for (j = 0; j < 4; j++)
1880 {
1881 real_t d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
1882 if (j == k) { d -= 1.0; }
1883 if (fabs(d) > 1.0e-12)
1884 {
1886 " k = " << k << ", j = " << j << ", d = " << d << endl;
1887 mfem_error();
1888 }
1889 }
1890 }
1891#endif
1892
1894 ip.x = ip.y = 0.0;
1895 Trans.SetIntPoint (&ip);
1896 // Trans must be linear (more to have embedding?)
1897 // set Jinv = |J| J^{-t} = adj(J)^t
1898 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
1899
1900 real_t vk[2];
1901 Vector xk (vk, 2);
1902
1903 for (k = 0; k < 4; k++)
1904 {
1905 Trans.Transform (Nodes.IntPoint (k), xk);
1906 ip.x = vk[0]; ip.y = vk[1];
1907 CalcVShape (ip, vshape);
1908 // vk = |J| J^{-t} nk
1909 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
1910 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
1911 for (j = 0; j < 4; j++)
1912 if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
1913 {
1914 I(k,j) = 0.0;
1915 }
1916 }
1917}
1918
1921 Vector &dofs) const
1922{
1923 real_t vk[2];
1924 Vector xk (vk, 2);
1925
1926 for (int k = 0; k < 4; k++)
1927 {
1928 Trans.SetIntPoint (&Nodes.IntPoint (k));
1929 // set Jinv = |J| J^{-t} = adj(J)^t
1930 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
1931
1932 vc.Eval (xk, Trans, Nodes.IntPoint (k));
1933 // xk^t |J| J^{-t} nk
1934 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
1935 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
1936 }
1937}
1938
1940 : VectorFiniteElement(2, Geometry::TRIANGLE, 8, 2, H_DIV)
1941{
1942 Nodes.IntPoint(0).x = 0.33333333333333333333;
1943 Nodes.IntPoint(0).y = 0.0;
1944 Nodes.IntPoint(1).x = 0.66666666666666666667;
1945 Nodes.IntPoint(1).y = 0.0;
1946 Nodes.IntPoint(2).x = 0.66666666666666666667;
1947 Nodes.IntPoint(2).y = 0.33333333333333333333;
1948 Nodes.IntPoint(3).x = 0.33333333333333333333;
1949 Nodes.IntPoint(3).y = 0.66666666666666666667;
1950 Nodes.IntPoint(4).x = 0.0;
1951 Nodes.IntPoint(4).y = 0.66666666666666666667;
1952 Nodes.IntPoint(5).x = 0.0;
1953 Nodes.IntPoint(5).y = 0.33333333333333333333;
1954 Nodes.IntPoint(6).x = 0.33333333333333333333;
1955 Nodes.IntPoint(6).y = 0.33333333333333333333;
1956 Nodes.IntPoint(7).x = 0.33333333333333333333;
1957 Nodes.IntPoint(7).y = 0.33333333333333333333;
1958}
1959
1961 DenseMatrix &shape) const
1962{
1963 real_t x = ip.x, y = ip.y;
1964
1965 shape(0,0) = -2 * x * (-1 + x + 2 * y);
1966 shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
1967 shape(1,0) = 2 * x * (x - y);
1968 shape(1,1) = 2 * (x - y) * (-1 + y);
1969 shape(2,0) = 2 * x * (-1 + 2 * x + y);
1970 shape(2,1) = 2 * y * (-1 + 2 * x + y);
1971 shape(3,0) = 2 * x * (-1 + x + 2 * y);
1972 shape(3,1) = 2 * y * (-1 + x + 2 * y);
1973 shape(4,0) = -2 * (-1 + x) * (x - y);
1974 shape(4,1) = 2 * y * (-x + y);
1975 shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
1976 shape(5,1) = -2 * y * (-1 + 2 * x + y);
1977 shape(6,0) = -3 * x * (-2 + 2 * x + y);
1978 shape(6,1) = -3 * y * (-1 + 2 * x + y);
1979 shape(7,0) = -3 * x * (-1 + x + 2 * y);
1980 shape(7,1) = -3 * y * (-2 + x + 2 * y);
1981}
1982
1984 Vector &divshape) const
1985{
1986 real_t x = ip.x, y = ip.y;
1987
1988 divshape(0) = -2 * (-4 + 3 * x + 6 * y);
1989 divshape(1) = 2 + 6 * x - 6 * y;
1990 divshape(2) = -4 + 12 * x + 6 * y;
1991 divshape(3) = -4 + 6 * x + 12 * y;
1992 divshape(4) = 2 - 6 * x + 6 * y;
1993 divshape(5) = -2 * (-4 + 6 * x + 3 * y);
1994 divshape(6) = -9 * (-1 + 2 * x + y);
1995 divshape(7) = -9 * (-1 + x + 2 * y);
1996}
1997
1998const real_t RT1TriangleFiniteElement::nk[8][2] =
1999{
2000 { 0,-1}, { 0,-1},
2001 { 1, 1}, { 1, 1},
2002 {-1, 0}, {-1, 0},
2003 { 1, 0}, { 0, 1}
2004};
2005
2007 ElementTransformation &Trans, DenseMatrix &I) const
2008{
2009 int k, j;
2012#endif
2013
2014#ifdef MFEM_DEBUG
2015 for (k = 0; k < 8; k++)
2016 {
2018 for (j = 0; j < 8; j++)
2019 {
2020 real_t d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
2021 if (j == k) { d -= 1.0; }
2022 if (fabs(d) > 1.0e-12)
2023 {
2025 " k = " << k << ", j = " << j << ", d = " << d << endl;
2026 mfem_error();
2027 }
2028 }
2029 }
2030#endif
2031
2033 ip.x = ip.y = 0.0;
2034 Trans.SetIntPoint (&ip);
2035 // Trans must be linear (more to have embedding?)
2036 // set Jinv = |J| J^{-t} = adj(J)^t
2037 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2038
2039 real_t vk[2];
2040 Vector xk (vk, 2);
2041
2042 for (k = 0; k < 8; k++)
2043 {
2044 Trans.Transform (Nodes.IntPoint (k), xk);
2045 ip.x = vk[0]; ip.y = vk[1];
2046 CalcVShape (ip, vshape);
2047 // vk = |J| J^{-t} nk
2048 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2049 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2050 for (j = 0; j < 8; j++)
2051 if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
2052 {
2053 I(k,j) = 0.0;
2054 }
2055 }
2056}
2057
2059 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
2060{
2061 real_t vk[2];
2062 Vector xk (vk, 2);
2063
2064 for (int k = 0; k < 8; k++)
2065 {
2066 Trans.SetIntPoint (&Nodes.IntPoint (k));
2067 // set Jinv = |J| J^{-t} = adj(J)^t
2068 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2069
2070 vc.Eval (xk, Trans, Nodes.IntPoint (k));
2071 // xk^t |J| J^{-t} nk
2072 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2073 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2074 dofs(k) *= 0.5;
2075 }
2076}
2077
2079 : VectorFiniteElement(2, Geometry::SQUARE, 12, 2, H_DIV,
2080 FunctionSpace::Qk)
2081{
2082 // y = 0
2083 Nodes.IntPoint(0).x = 1./3.;
2084 Nodes.IntPoint(0).y = 0.0;
2085 Nodes.IntPoint(1).x = 2./3.;
2086 Nodes.IntPoint(1).y = 0.0;
2087 // x = 1
2088 Nodes.IntPoint(2).x = 1.0;
2089 Nodes.IntPoint(2).y = 1./3.;
2090 Nodes.IntPoint(3).x = 1.0;
2091 Nodes.IntPoint(3).y = 2./3.;
2092 // y = 1
2093 Nodes.IntPoint(4).x = 2./3.;
2094 Nodes.IntPoint(4).y = 1.0;
2095 Nodes.IntPoint(5).x = 1./3.;
2096 Nodes.IntPoint(5).y = 1.0;
2097 // x = 0
2098 Nodes.IntPoint(6).x = 0.0;
2099 Nodes.IntPoint(6).y = 2./3.;
2100 Nodes.IntPoint(7).x = 0.0;
2101 Nodes.IntPoint(7).y = 1./3.;
2102 // x = 0.5 (interior)
2103 Nodes.IntPoint(8).x = 0.5;
2104 Nodes.IntPoint(8).y = 1./3.;
2105 Nodes.IntPoint(9).x = 0.5;
2106 Nodes.IntPoint(9).y = 2./3.;
2107 // y = 0.5 (interior)
2108 Nodes.IntPoint(10).x = 1./3.;
2109 Nodes.IntPoint(10).y = 0.5;
2110 Nodes.IntPoint(11).x = 2./3.;
2111 Nodes.IntPoint(11).y = 0.5;
2112}
2113
2115 DenseMatrix &shape) const
2116{
2117 real_t x = ip.x, y = ip.y;
2118
2119 // y = 0
2120 shape(0,0) = 0;
2121 shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
2122 shape(1,0) = 0;
2123 shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
2124 // x = 1
2125 shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
2126 shape(2,1) = 0;
2127 shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
2128 shape(3,1) = 0;
2129 // y = 1
2130 shape(4,0) = 0;
2131 shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
2132 shape(5,0) = 0;
2133 shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
2134 // x = 0
2135 shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
2136 shape(6,1) = 0;
2137 shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
2138 shape(7,1) = 0;
2139 // x = 0.5 (interior)
2140 shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
2141 shape(8,1) = 0;
2142 shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
2143 shape(9,1) = 0;
2144 // y = 0.5 (interior)
2145 shape(10,0) = 0;
2146 shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
2147 shape(11,0) = 0;
2148 shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
2149}
2150
2152 Vector &divshape) const
2153{
2154 real_t x = ip.x, y = ip.y;
2155
2156 divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
2157 divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
2158 divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
2159 divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
2160 divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
2161 divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
2162 divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
2163 divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
2164 divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
2165 divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
2166 divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
2167 divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
2168}
2169
2171{
2172 // y = 0
2173 {0,-1}, {0,-1},
2174 // X = 1
2175 {1, 0}, {1, 0},
2176 // y = 1
2177 {0, 1}, {0, 1},
2178 // x = 0
2179 {-1,0}, {-1,0},
2180 // x = 0.5 (interior)
2181 {1, 0}, {1, 0},
2182 // y = 0.5 (interior)
2183 {0, 1}, {0, 1}
2184};
2185
2187 ElementTransformation &Trans, DenseMatrix &I) const
2188{
2189 int k, j;
2192#endif
2193
2194#ifdef MFEM_DEBUG
2195 for (k = 0; k < 12; k++)
2196 {
2198 for (j = 0; j < 12; j++)
2199 {
2200 real_t d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
2201 if (j == k) { d -= 1.0; }
2202 if (fabs(d) > 1.0e-12)
2203 {
2205 " k = " << k << ", j = " << j << ", d = " << d << endl;
2206 mfem_error();
2207 }
2208 }
2209 }
2210#endif
2211
2213 ip.x = ip.y = 0.0;
2214 Trans.SetIntPoint (&ip);
2215 // Trans must be linear (more to have embedding?)
2216 // set Jinv = |J| J^{-t} = adj(J)^t
2217 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2218
2219 real_t vk[2];
2220 Vector xk (vk, 2);
2221
2222 for (k = 0; k < 12; k++)
2223 {
2224 Trans.Transform (Nodes.IntPoint (k), xk);
2225 ip.x = vk[0]; ip.y = vk[1];
2226 CalcVShape (ip, vshape);
2227 // vk = |J| J^{-t} nk
2228 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2229 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2230 for (j = 0; j < 12; j++)
2231 if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
2232 {
2233 I(k,j) = 0.0;
2234 }
2235 }
2236}
2237
2239 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
2240{
2241 real_t vk[2];
2242 Vector xk (vk, 2);
2243
2244 for (int k = 0; k < 12; k++)
2245 {
2246 Trans.SetIntPoint (&Nodes.IntPoint (k));
2247 // set Jinv = |J| J^{-t} = adj(J)^t
2248 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2249
2250 vc.Eval (xk, Trans, Nodes.IntPoint (k));
2251 // xk^t |J| J^{-t} nk
2252 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2253 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2254 }
2255}
2256
2257const real_t RT2TriangleFiniteElement::M[15][15] =
2258{
2259 // *INDENT-OFF*
2260 {
2261 0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
2262 0, 24.442740046346700787, -16.647580015448900262, -12.,
2263 -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
2264 12., 30.590320061795601049, 15.295160030897800524
2265 },
2266 {
2267 0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
2268 -15., 10.5
2269 },
2270 {
2271 0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
2272 0, -3.4427400463467007866, -7.3524199845510997378, -12.,
2273 4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
2274 12., -6.5903200617956010489, -3.2951600308978005244
2275 },
2276 {
2277 0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
2278 -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
2279 0, -3.2951600308978005244
2280 },
2281 {
2282 0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
2283 36., 10.5
2284 },
2285 {
2286 0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
2287 2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
2288 0, 15.295160030897800524
2289 },
2290 {
2291 -0.67620999227554986889, 0, -3.4427400463467007866, 0,
2292 7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
2293 -0.76209992275549868892, 4.1189500386222506555, -12.,
2294 -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
2295 12.
2296 },
2297 {
2298 1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
2299 -15., -15.
2300 },
2301 {
2302 -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
2303 5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
2304 -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
2305 30.590320061795601049, 12.
2306 },
2307 { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
2308 { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
2309 { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
2310 { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
2311 { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
2312 { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
2313 // *INDENT-ON*
2314};
2315
2317 : VectorFiniteElement(2, Geometry::TRIANGLE, 15, 3, H_DIV)
2318{
2319 const real_t p = 0.11270166537925831148;
2320
2321 Nodes.IntPoint(0).x = p;
2322 Nodes.IntPoint(0).y = 0.0;
2323 Nodes.IntPoint(1).x = 0.5;
2324 Nodes.IntPoint(1).y = 0.0;
2325 Nodes.IntPoint(2).x = 1.-p;
2326 Nodes.IntPoint(2).y = 0.0;
2327 Nodes.IntPoint(3).x = 1.-p;
2328 Nodes.IntPoint(3).y = p;
2329 Nodes.IntPoint(4).x = 0.5;
2330 Nodes.IntPoint(4).y = 0.5;
2331 Nodes.IntPoint(5).x = p;
2332 Nodes.IntPoint(5).y = 1.-p;
2333 Nodes.IntPoint(6).x = 0.0;
2334 Nodes.IntPoint(6).y = 1.-p;
2335 Nodes.IntPoint(7).x = 0.0;
2336 Nodes.IntPoint(7).y = 0.5;
2337 Nodes.IntPoint(8).x = 0.0;
2338 Nodes.IntPoint(8).y = p;
2339 Nodes.IntPoint(9).x = 0.25;
2340 Nodes.IntPoint(9).y = 0.25;
2341 Nodes.IntPoint(10).x = 0.25;
2342 Nodes.IntPoint(10).y = 0.25;
2343 Nodes.IntPoint(11).x = 0.5;
2344 Nodes.IntPoint(11).y = 0.25;
2345 Nodes.IntPoint(12).x = 0.5;
2346 Nodes.IntPoint(12).y = 0.25;
2347 Nodes.IntPoint(13).x = 0.25;
2348 Nodes.IntPoint(13).y = 0.5;
2349 Nodes.IntPoint(14).x = 0.25;
2350 Nodes.IntPoint(14).y = 0.5;
2351}
2352
2354 DenseMatrix &shape) const
2355{
2356 real_t x = ip.x, y = ip.y;
2357
2358 real_t Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
2359 x*x*y, x*y*y
2360 };
2361 real_t By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
2362 x*x*y, x*y*y, y*y*y
2363 };
2364
2365 for (int i = 0; i < 15; i++)
2366 {
2367 real_t cx = 0.0, cy = 0.0;
2368 for (int j = 0; j < 15; j++)
2369 {
2370 cx += M[i][j] * Bx[j];
2371 cy += M[i][j] * By[j];
2372 }
2373 shape(i,0) = cx;
2374 shape(i,1) = cy;
2375 }
2376}
2377
2379 Vector &divshape) const
2380{
2381 real_t x = ip.x, y = ip.y;
2382 constexpr real_t f2 = 2.0;
2383 constexpr real_t f4 = 4.0;
2384
2385 real_t DivB[15] = {0., 0., 1., 0., 0., 1., f2*x, 0., y, x, 0., f2*y,
2386 f4*x*x, f4*x*y, f4*y*y
2387 };
2388
2389 for (int i = 0; i < 15; i++)
2390 {
2391 real_t div = 0.0;
2392 for (int j = 0; j < 15; j++)
2393 {
2394 div += M[i][j] * DivB[j];
2395 }
2396 divshape(i) = div;
2397 }
2398}
2399
2401
2403
2405 : VectorFiniteElement(2, Geometry::SQUARE, 24, 3, H_DIV,
2406 FunctionSpace::Qk)
2407{
2408 // y = 0 (pt[0])
2409 Nodes.IntPoint(0).x = dpt[0]; Nodes.IntPoint(0).y = pt[0];
2410 Nodes.IntPoint(1).x = dpt[1]; Nodes.IntPoint(1).y = pt[0];
2411 Nodes.IntPoint(2).x = dpt[2]; Nodes.IntPoint(2).y = pt[0];
2412 // x = 1 (pt[3])
2413 Nodes.IntPoint(3).x = pt[3]; Nodes.IntPoint(3).y = dpt[0];
2414 Nodes.IntPoint(4).x = pt[3]; Nodes.IntPoint(4).y = dpt[1];
2415 Nodes.IntPoint(5).x = pt[3]; Nodes.IntPoint(5).y = dpt[2];
2416 // y = 1 (pt[3])
2417 Nodes.IntPoint(6).x = dpt[2]; Nodes.IntPoint(6).y = pt[3];
2418 Nodes.IntPoint(7).x = dpt[1]; Nodes.IntPoint(7).y = pt[3];
2419 Nodes.IntPoint(8).x = dpt[0]; Nodes.IntPoint(8).y = pt[3];
2420 // x = 0 (pt[0])
2421 Nodes.IntPoint(9).x = pt[0]; Nodes.IntPoint(9).y = dpt[2];
2422 Nodes.IntPoint(10).x = pt[0]; Nodes.IntPoint(10).y = dpt[1];
2423 Nodes.IntPoint(11).x = pt[0]; Nodes.IntPoint(11).y = dpt[0];
2424 // x = pt[1] (interior)
2425 Nodes.IntPoint(12).x = pt[1]; Nodes.IntPoint(12).y = dpt[0];
2426 Nodes.IntPoint(13).x = pt[1]; Nodes.IntPoint(13).y = dpt[1];
2427 Nodes.IntPoint(14).x = pt[1]; Nodes.IntPoint(14).y = dpt[2];
2428 // x = pt[2] (interior)
2429 Nodes.IntPoint(15).x = pt[2]; Nodes.IntPoint(15).y = dpt[0];
2430 Nodes.IntPoint(16).x = pt[2]; Nodes.IntPoint(16).y = dpt[1];
2431 Nodes.IntPoint(17).x = pt[2]; Nodes.IntPoint(17).y = dpt[2];
2432 // y = pt[1] (interior)
2433 Nodes.IntPoint(18).x = dpt[0]; Nodes.IntPoint(18).y = pt[1];
2434 Nodes.IntPoint(19).x = dpt[1]; Nodes.IntPoint(19).y = pt[1];
2435 Nodes.IntPoint(20).x = dpt[2]; Nodes.IntPoint(20).y = pt[1];
2436 // y = pt[2] (interior)
2437 Nodes.IntPoint(21).x = dpt[0]; Nodes.IntPoint(21).y = pt[2];
2438 Nodes.IntPoint(22).x = dpt[1]; Nodes.IntPoint(22).y = pt[2];
2439 Nodes.IntPoint(23).x = dpt[2]; Nodes.IntPoint(23).y = pt[2];
2440}
2441
2443 DenseMatrix &shape) const
2444{
2445 real_t x = ip.x, y = ip.y;
2446
2447 real_t ax0 = pt[0] - x;
2448 real_t ax1 = pt[1] - x;
2449 real_t ax2 = pt[2] - x;
2450 real_t ax3 = pt[3] - x;
2451
2452 real_t by0 = dpt[0] - y;
2453 real_t by1 = dpt[1] - y;
2454 real_t by2 = dpt[2] - y;
2455
2456 real_t ay0 = pt[0] - y;
2457 real_t ay1 = pt[1] - y;
2458 real_t ay2 = pt[2] - y;
2459 real_t ay3 = pt[3] - y;
2460
2461 real_t bx0 = dpt[0] - x;
2462 real_t bx1 = dpt[1] - x;
2463 real_t bx2 = dpt[2] - x;
2464
2465 real_t A01 = pt[0] - pt[1];
2466 real_t A02 = pt[0] - pt[2];
2467 real_t A12 = pt[1] - pt[2];
2468 real_t A03 = pt[0] - pt[3];
2469 real_t A13 = pt[1] - pt[3];
2470 real_t A23 = pt[2] - pt[3];
2471
2472 real_t B01 = dpt[0] - dpt[1];
2473 real_t B02 = dpt[0] - dpt[2];
2474 real_t B12 = dpt[1] - dpt[2];
2475
2476 real_t tx0 = (bx1*bx2)/(B01*B02);
2477 real_t tx1 = -(bx0*bx2)/(B01*B12);
2478 real_t tx2 = (bx0*bx1)/(B02*B12);
2479
2480 real_t ty0 = (by1*by2)/(B01*B02);
2481 real_t ty1 = -(by0*by2)/(B01*B12);
2482 real_t ty2 = (by0*by1)/(B02*B12);
2483
2484 // y = 0 (p[0])
2485 shape(0, 0) = 0;
2486 shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
2487 shape(1, 0) = 0;
2488 shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
2489 shape(2, 0) = 0;
2490 shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
2491 // x = 1 (p[3])
2492 shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
2493 shape(3, 1) = 0;
2494 shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
2495 shape(4, 1) = 0;
2496 shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
2497 shape(5, 1) = 0;
2498 // y = 1 (p[3])
2499 shape(6, 0) = 0;
2500 shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
2501 shape(7, 0) = 0;
2502 shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
2503 shape(8, 0) = 0;
2504 shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
2505 // x = 0 (p[0])
2506 shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
2507 shape(9, 1) = 0;
2508 shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
2509 shape(10, 1) = 0;
2510 shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
2511 shape(11, 1) = 0;
2512 // x = p[1] (interior)
2513 shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
2514 shape(12, 1) = 0;
2515 shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
2516 shape(13, 1) = 0;
2517 shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
2518 shape(14, 1) = 0;
2519 // x = p[2] (interior)
2520 shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
2521 shape(15, 1) = 0;
2522 shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
2523 shape(16, 1) = 0;
2524 shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
2525 shape(17, 1) = 0;
2526 // y = p[1] (interior)
2527 shape(18, 0) = 0;
2528 shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
2529 shape(19, 0) = 0;
2530 shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
2531 shape(20, 0) = 0;
2532 shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
2533 // y = p[2] (interior)
2534 shape(21, 0) = 0;
2535 shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
2536 shape(22, 0) = 0;
2537 shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
2538 shape(23, 0) = 0;
2539 shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
2540}
2541
2543 Vector &divshape) const
2544{
2545 real_t x = ip.x, y = ip.y;
2546
2547 real_t a01 = pt[0]*pt[1];
2548 real_t a02 = pt[0]*pt[2];
2549 real_t a12 = pt[1]*pt[2];
2550 real_t a03 = pt[0]*pt[3];
2551 real_t a13 = pt[1]*pt[3];
2552 real_t a23 = pt[2]*pt[3];
2553
2554 real_t bx0 = dpt[0] - x;
2555 real_t bx1 = dpt[1] - x;
2556 real_t bx2 = dpt[2] - x;
2557
2558 real_t by0 = dpt[0] - y;
2559 real_t by1 = dpt[1] - y;
2560 real_t by2 = dpt[2] - y;
2561
2562 real_t A01 = pt[0] - pt[1];
2563 real_t A02 = pt[0] - pt[2];
2564 real_t A12 = pt[1] - pt[2];
2565 real_t A03 = pt[0] - pt[3];
2566 real_t A13 = pt[1] - pt[3];
2567 real_t A23 = pt[2] - pt[3];
2568
2569 real_t A012 = pt[0] + pt[1] + pt[2];
2570 real_t A013 = pt[0] + pt[1] + pt[3];
2571 real_t A023 = pt[0] + pt[2] + pt[3];
2572 real_t A123 = pt[1] + pt[2] + pt[3];
2573
2574 real_t B01 = dpt[0] - dpt[1];
2575 real_t B02 = dpt[0] - dpt[2];
2576 real_t B12 = dpt[1] - dpt[2];
2577
2578 real_t tx0 = (bx1*bx2)/(B01*B02);
2579 real_t tx1 = -(bx0*bx2)/(B01*B12);
2580 real_t tx2 = (bx0*bx1)/(B02*B12);
2581
2582 real_t ty0 = (by1*by2)/(B01*B02);
2583 real_t ty1 = -(by0*by2)/(B01*B12);
2584 real_t ty2 = (by0*by1)/(B02*B12);
2585
2586 // y = 0 (p[0])
2587 divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
2588 divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
2589 divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
2590 // x = 1 (p[3])
2591 divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
2592 divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
2593 divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
2594 // y = 1 (p[3])
2595 divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
2596 divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
2597 divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
2598 // x = 0 (p[0])
2599 divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
2600 divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
2601 divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
2602 // x = p[1] (interior)
2603 divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
2604 divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
2605 divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
2606 // x = p[2] (interior)
2607 divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
2608 divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
2609 divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
2610 // y = p[1] (interior)
2611 divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
2612 divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
2613 divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
2614 // y = p[2] (interior)
2615 divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
2616 divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
2617 divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
2618}
2619
2621{
2622 // y = 0
2623 {0,-1}, {0,-1}, {0,-1},
2624 // x = 1
2625 {1, 0}, {1, 0}, {1, 0},
2626 // y = 1
2627 {0, 1}, {0, 1}, {0, 1},
2628 // x = 0
2629 {-1,0}, {-1,0}, {-1,0},
2630 // x = p[1] (interior)
2631 {1, 0}, {1, 0}, {1, 0},
2632 // x = p[2] (interior)
2633 {1, 0}, {1, 0}, {1, 0},
2634 // y = p[1] (interior)
2635 {0, 1}, {0, 1}, {0, 1},
2636 // y = p[1] (interior)
2637 {0, 1}, {0, 1}, {0, 1}
2638};
2639
2641 ElementTransformation &Trans, DenseMatrix &I) const
2642{
2643 int k, j;
2646#endif
2647
2648#ifdef MFEM_DEBUG
2649 for (k = 0; k < 24; k++)
2650 {
2652 for (j = 0; j < 24; j++)
2653 {
2654 real_t d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
2655 if (j == k) { d -= 1.0; }
2656 if (fabs(d) > 1.0e-12)
2657 {
2659 " k = " << k << ", j = " << j << ", d = " << d << endl;
2660 mfem_error();
2661 }
2662 }
2663 }
2664#endif
2665
2667 ip.x = ip.y = 0.0;
2668 Trans.SetIntPoint (&ip);
2669 // Trans must be linear (more to have embedding?)
2670 // set Jinv = |J| J^{-t} = adj(J)^t
2671 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2672
2673 real_t vk[2];
2674 Vector xk (vk, 2);
2675
2676 for (k = 0; k < 24; k++)
2677 {
2678 Trans.Transform (Nodes.IntPoint (k), xk);
2679 ip.x = vk[0]; ip.y = vk[1];
2680 CalcVShape (ip, vshape);
2681 // vk = |J| J^{-t} nk
2682 vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2683 vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2684 for (j = 0; j < 24; j++)
2685 if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
2686 {
2687 I(k,j) = 0.0;
2688 }
2689 }
2690}
2691
2693 VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
2694{
2695 real_t vk[2];
2696 Vector xk (vk, 2);
2697
2698 for (int k = 0; k < 24; k++)
2699 {
2700 Trans.SetIntPoint (&Nodes.IntPoint (k));
2701 // set Jinv = |J| J^{-t} = adj(J)^t
2702 const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2703
2704 vc.Eval (xk, Trans, Nodes.IntPoint (k));
2705 // xk^t |J| J^{-t} nk
2706 dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2707 vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2708 }
2709}
2710
2712 : NodalFiniteElement(1, Geometry::SEGMENT, 2, 1)
2713{
2714 Nodes.IntPoint(0).x = 0.33333333333333333333;
2715 Nodes.IntPoint(1).x = 0.66666666666666666667;
2716}
2717
2719 Vector &shape) const
2720{
2721 real_t x = ip.x;
2722
2723 shape(0) = 2. - 3. * x;
2724 shape(1) = 3. * x - 1.;
2725}
2726
2728 DenseMatrix &dshape) const
2729{
2730 dshape(0,0) = -3.;
2731 dshape(1,0) = 3.;
2732}
2733
2734
2736 : NodalFiniteElement(1, Geometry::SEGMENT, 3, 2)
2737{
2738 const real_t p = 0.11270166537925831148;
2739
2740 Nodes.IntPoint(0).x = p;
2741 Nodes.IntPoint(1).x = 0.5;
2742 Nodes.IntPoint(2).x = 1.-p;
2743}
2744
2746 Vector &shape) const
2747{
2748 const real_t p = 0.11270166537925831148;
2749 const real_t w = 1./((1-2*p)*(1-2*p));
2750 real_t x = ip.x;
2751
2752 shape(0) = (2*x-1)*(x-1+p)*w;
2753 shape(1) = 4*(x-1+p)*(p-x)*w;
2754 shape(2) = (2*x-1)*(x-p)*w;
2755}
2756
2758 DenseMatrix &dshape) const
2759{
2760 const real_t p = 0.11270166537925831148;
2761 const real_t w = 1./((1-2*p)*(1-2*p));
2762 real_t x = ip.x;
2763
2764 dshape(0,0) = (-3+4*x+2*p)*w;
2765 dshape(1,0) = (4-8*x)*w;
2766 dshape(2,0) = (-1+4*x-2*p)*w;
2767}
2768
2769
2771 : NodalFiniteElement(1, Geometry::SEGMENT, degree+1, degree)
2772{
2773 int i, m = degree;
2774
2775 Nodes.IntPoint(0).x = 0.0;
2776 Nodes.IntPoint(1).x = 1.0;
2777 for (i = 1; i < m; i++)
2778 {
2779 Nodes.IntPoint(i+1).x = real_t(i) / m;
2780 }
2781
2782 rwk.SetSize(degree+1);
2784 rxxk.SetSize(degree+1);
2785#endif
2786
2787 rwk(0) = 1.0;
2788 for (i = 1; i <= m; i++)
2789 {
2790 rwk(i) = rwk(i-1) * ( (real_t)(m) / (real_t)(i) );
2791 }
2792 for (i = 0; i < m/2+1; i++)
2793 {
2794 rwk(m-i) = ( rwk(i) *= rwk(m-i) );
2795 }
2796 for (i = m-1; i >= 0; i -= 2)
2797 {
2798 rwk(i) = -rwk(i);
2799 }
2800}
2801
2803 Vector &shape) const
2804{
2805 real_t w, wk, x = ip.x;
2806 int i, k, m = GetOrder();
2807
2809 Vector rxxk(m+1);
2810#endif
2811
2812 k = (int) floor ( m * x + 0.5 );
2813 k = k > m ? m : k < 0 ? 0 : k; // clamp k to [0,m]
2814
2815 wk = 1.0;
2816 for (i = 0; i <= m; i++)
2817 if (i != k)
2818 {
2819 wk *= ( rxxk(i) = x - (real_t)(i) / m );
2820 }
2821 w = wk * ( rxxk(k) = x - (real_t)(k) / m );
2822
2823 if (k != 0)
2824 {
2825 shape(0) = w * rwk(0) / rxxk(0);
2826 }
2827 else
2828 {
2829 shape(0) = wk * rwk(0);
2830 }
2831 if (k != m)
2832 {
2833 shape(1) = w * rwk(m) / rxxk(m);
2834 }
2835 else
2836 {
2837 shape(1) = wk * rwk(k);
2838 }
2839 for (i = 1; i < m; i++)
2840 if (i != k)
2841 {
2842 shape(i+1) = w * rwk(i) / rxxk(i);
2843 }
2844 else
2845 {
2846 shape(k+1) = wk * rwk(k);
2847 }
2848}
2849
2851 DenseMatrix &dshape) const
2852{
2853 real_t s, srx, w, wk, x = ip.x;
2854 int i, k, m = GetOrder();
2855
2857 Vector rxxk(m+1);
2858#endif
2859
2860 k = (int) floor ( m * x + 0.5 );
2861 k = k > m ? m : k < 0 ? 0 : k; // clamp k to [0,m]
2862
2863 wk = 1.0;
2864 for (i = 0; i <= m; i++)
2865 if (i != k)
2866 {
2867 wk *= ( rxxk(i) = x - (real_t)(i) / m );
2868 }
2869 w = wk * ( rxxk(k) = x - (real_t)(k) / m );
2870
2871 for (i = 0; i <= m; i++)
2872 {
2873 rxxk(i) = 1.0 / rxxk(i);
2874 }
2875 srx = 0.0;
2876 for (i = 0; i <= m; i++)
2877 if (i != k)
2878 {
2879 srx += rxxk(i);
2880 }
2881 s = w * srx + wk;
2882
2883 if (k != 0)
2884 {
2885 dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
2886 }
2887 else
2888 {
2889 dshape(0,0) = wk * srx * rwk(0);
2890 }
2891 if (k != m)
2892 {
2893 dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
2894 }
2895 else
2896 {
2897 dshape(1,0) = wk * srx * rwk(k);
2898 }
2899 for (i = 1; i < m; i++)
2900 if (i != k)
2901 {
2902 dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
2903 }
2904 else
2905 {
2906 dshape(k+1,0) = wk * srx * rwk(k);
2907 }
2908}
2909
2910
2912 : NodalFiniteElement(3, Geometry::TETRAHEDRON, 4, 1)
2913{
2914 Nodes.IntPoint(0).x = 0.33333333333333333333;
2915 Nodes.IntPoint(0).y = 0.33333333333333333333;
2916 Nodes.IntPoint(0).z = 0.33333333333333333333;
2917
2918 Nodes.IntPoint(1).x = 0.0;
2919 Nodes.IntPoint(1).y = 0.33333333333333333333;
2920 Nodes.IntPoint(1).z = 0.33333333333333333333;
2921
2922 Nodes.IntPoint(2).x = 0.33333333333333333333;
2923 Nodes.IntPoint(2).y = 0.0;
2924 Nodes.IntPoint(2).z = 0.33333333333333333333;
2925
2926 Nodes.IntPoint(3).x = 0.33333333333333333333;
2927 Nodes.IntPoint(3).y = 0.33333333333333333333;
2928 Nodes.IntPoint(3).z = 0.0;
2929
2930}
2931
2933 Vector &shape) const
2934{
2935 real_t L0, L1, L2, L3;
2936
2937 L1 = ip.x; L2 = ip.y; L3 = ip.z; L0 = 1.0 - L1 - L2 - L3;
2938 shape(0) = 1.0 - 3.0 * L0;
2939 shape(1) = 1.0 - 3.0 * L1;
2940 shape(2) = 1.0 - 3.0 * L2;
2941 shape(3) = 1.0 - 3.0 * L3;
2942}
2943
2945 DenseMatrix &dshape) const
2946{
2947 dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
2948 dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2949 dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
2950 dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
2951}
2952
2953
2955 : NodalFiniteElement(3, Geometry::TETRAHEDRON, 1, 0)
2956{
2957 Nodes.IntPoint(0).x = 0.25;
2958 Nodes.IntPoint(0).y = 0.25;
2959 Nodes.IntPoint(0).z = 0.25;
2960}
2961
2963 Vector &shape) const
2964{
2965 shape(0) = 1.0;
2966}
2967
2969 DenseMatrix &dshape) const
2970{
2971 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
2972}
2973
2974
2976 : NodalFiniteElement(3, Geometry::CUBE, 1, 0, FunctionSpace::Qk)
2977{
2978 Nodes.IntPoint(0).x = 0.5;
2979 Nodes.IntPoint(0).y = 0.5;
2980 Nodes.IntPoint(0).z = 0.5;
2981}
2982
2984 Vector &shape) const
2985{
2986 shape(0) = 1.0;
2987}
2988
2990 DenseMatrix &dshape) const
2991{
2992 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
2993}
2994
2995
2997 : NodalFiniteElement(3, Geometry::PRISM, 1, 0, FunctionSpace::Qk)
2998{
2999 Nodes.IntPoint(0).x = 1.0 / 3.0;
3000 Nodes.IntPoint(0).y = 1.0 / 3.0;
3001 Nodes.IntPoint(0).z = 0.5;
3002}
3003
3005 Vector &shape) const
3006{
3007 shape(0) = 1.0;
3008}
3009
3011 DenseMatrix &dshape) const
3012{
3013 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3014}
3015
3016
3018 : NodalFiniteElement(3, Geometry::PYRAMID, 1, 0, FunctionSpace::Qk)
3019{
3020 Nodes.IntPoint(0).x = 0.375;
3021 Nodes.IntPoint(0).y = 0.375;
3022 Nodes.IntPoint(0).z = 0.25;
3023}
3024
3026 Vector &shape) const
3027{
3028 shape(0) = 1.0;
3029}
3030
3032 DenseMatrix &dshape) const
3033{
3034 dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3035}
3036
3037
3039 : NodalFiniteElement(3, Geometry::CUBE, (degree+1)*(degree+1)*(degree+1),
3040 degree, FunctionSpace::Qk)
3041{
3042 if (degree == 2)
3043 {
3044 I = new int[dof];
3045 J = new int[dof];
3046 K = new int[dof];
3047 // nodes
3048 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3049 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3050 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3051 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3052 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3053 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3054 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3055 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3056 // edges
3057 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3058 I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3059 I[10] = 2; J[10] = 1; K[10] = 0;
3060 I[11] = 0; J[11] = 2; K[11] = 0;
3061 I[12] = 2; J[12] = 0; K[12] = 1;
3062 I[13] = 1; J[13] = 2; K[13] = 1;
3063 I[14] = 2; J[14] = 1; K[14] = 1;
3064 I[15] = 0; J[15] = 2; K[15] = 1;
3065 I[16] = 0; J[16] = 0; K[16] = 2;
3066 I[17] = 1; J[17] = 0; K[17] = 2;
3067 I[18] = 1; J[18] = 1; K[18] = 2;
3068 I[19] = 0; J[19] = 1; K[19] = 2;
3069 // faces
3070 I[20] = 2; J[20] = 2; K[20] = 0;
3071 I[21] = 2; J[21] = 0; K[21] = 2;
3072 I[22] = 1; J[22] = 2; K[22] = 2;
3073 I[23] = 2; J[23] = 1; K[23] = 2;
3074 I[24] = 0; J[24] = 2; K[24] = 2;
3075 I[25] = 2; J[25] = 2; K[25] = 1;
3076 // element
3077 I[26] = 2; J[26] = 2; K[26] = 2;
3078 }
3079 else if (degree == 3)
3080 {
3081 I = new int[dof];
3082 J = new int[dof];
3083 K = new int[dof];
3084 // nodes
3085 I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3086 I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3087 I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3088 I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3089 I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3090 I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3091 I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3092 I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3093 // edges
3094 I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3095 I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3096 I[10] = 1; J[10] = 2; K[10] = 0;
3097 I[11] = 1; J[11] = 3; K[11] = 0;
3098 I[12] = 2; J[12] = 1; K[12] = 0;
3099 I[13] = 3; J[13] = 1; K[13] = 0;
3100 I[14] = 0; J[14] = 2; K[14] = 0;
3101 I[15] = 0; J[15] = 3; K[15] = 0;
3102 I[16] = 2; J[16] = 0; K[16] = 1;
3103 I[17] = 3; J[17] = 0; K[17] = 1;
3104 I[18] = 1; J[18] = 2; K[18] = 1;
3105 I[19] = 1; J[19] = 3; K[19] = 1;
3106 I[20] = 2; J[20] = 1; K[20] = 1;
3107 I[21] = 3; J[21] = 1; K[21] = 1;
3108 I[22] = 0; J[22] = 2; K[22] = 1;
3109 I[23] = 0; J[23] = 3; K[23] = 1;
3110 I[24] = 0; J[24] = 0; K[24] = 2;
3111 I[25] = 0; J[25] = 0; K[25] = 3;
3112 I[26] = 1; J[26] = 0; K[26] = 2;
3113 I[27] = 1; J[27] = 0; K[27] = 3;
3114 I[28] = 1; J[28] = 1; K[28] = 2;
3115 I[29] = 1; J[29] = 1; K[29] = 3;
3116 I[30] = 0; J[30] = 1; K[30] = 2;
3117 I[31] = 0; J[31] = 1; K[31] = 3;
3118 // faces
3119 I[32] = 2; J[32] = 3; K[32] = 0;
3120 I[33] = 3; J[33] = 3; K[33] = 0;
3121 I[34] = 2; J[34] = 2; K[34] = 0;
3122 I[35] = 3; J[35] = 2; K[35] = 0;
3123 I[36] = 2; J[36] = 0; K[36] = 2;
3124 I[37] = 3; J[37] = 0; K[37] = 2;
3125 I[38] = 2; J[38] = 0; K[38] = 3;
3126 I[39] = 3; J[39] = 0; K[39] = 3;
3127 I[40] = 1; J[40] = 2; K[40] = 2;
3128 I[41] = 1; J[41] = 3; K[41] = 2;
3129 I[42] = 1; J[42] = 2; K[42] = 3;
3130 I[43] = 1; J[43] = 3; K[43] = 3;
3131 I[44] = 3; J[44] = 1; K[44] = 2;
3132 I[45] = 2; J[45] = 1; K[45] = 2;
3133 I[46] = 3; J[46] = 1; K[46] = 3;
3134 I[47] = 2; J[47] = 1; K[47] = 3;
3135 I[48] = 0; J[48] = 3; K[48] = 2;
3136 I[49] = 0; J[49] = 2; K[49] = 2;
3137 I[50] = 0; J[50] = 3; K[50] = 3;
3138 I[51] = 0; J[51] = 2; K[51] = 3;
3139 I[52] = 2; J[52] = 2; K[52] = 1;
3140 I[53] = 3; J[53] = 2; K[53] = 1;
3141 I[54] = 2; J[54] = 3; K[54] = 1;
3142 I[55] = 3; J[55] = 3; K[55] = 1;
3143 // element
3144 I[56] = 2; J[56] = 2; K[56] = 2;
3145 I[57] = 3; J[57] = 2; K[57] = 2;
3146 I[58] = 3; J[58] = 3; K[58] = 2;
3147 I[59] = 2; J[59] = 3; K[59] = 2;
3148 I[60] = 2; J[60] = 2; K[60] = 3;
3149 I[61] = 3; J[61] = 2; K[61] = 3;
3150 I[62] = 3; J[62] = 3; K[62] = 3;
3151 I[63] = 2; J[63] = 3; K[63] = 3;
3152 }
3153 else
3154 {
3155 mfem_error ("LagrangeHexFiniteElement::LagrangeHexFiniteElement");
3156 }
3157
3158 fe1d = new Lagrange1DFiniteElement(degree);
3159 dof1d = fe1d -> GetDof();
3160
3162 shape1dx.SetSize(dof1d);
3163 shape1dy.SetSize(dof1d);
3164 shape1dz.SetSize(dof1d);
3165
3166 dshape1dx.SetSize(dof1d,1);
3167 dshape1dy.SetSize(dof1d,1);
3168 dshape1dz.SetSize(dof1d,1);
3169#endif
3170
3171 for (int n = 0; n < dof; n++)
3172 {
3173 Nodes.IntPoint(n).x = fe1d -> GetNodes().IntPoint(I[n]).x;
3174 Nodes.IntPoint(n).y = fe1d -> GetNodes().IntPoint(J[n]).x;
3175 Nodes.IntPoint(n).z = fe1d -> GetNodes().IntPoint(K[n]).x;
3176 }
3177}
3178
3180 Vector &shape) const
3181{
3182 IntegrationPoint ipy, ipz;
3183 ipy.x = ip.y;
3184 ipz.x = ip.z;
3185
3187 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3188#endif
3189
3190 fe1d -> CalcShape(ip, shape1dx);
3191 fe1d -> CalcShape(ipy, shape1dy);
3192 fe1d -> CalcShape(ipz, shape1dz);
3193
3194 for (int n = 0; n < dof; n++)
3195 {
3196 shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
3197 }
3198}
3199
3201 DenseMatrix &dshape) const
3202{
3203 IntegrationPoint ipy, ipz;
3204 ipy.x = ip.y;
3205 ipz.x = ip.z;
3206
3208 Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3209 DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
3210#endif
3211
3212 fe1d -> CalcShape(ip, shape1dx);
3213 fe1d -> CalcShape(ipy, shape1dy);
3214 fe1d -> CalcShape(ipz, shape1dz);
3215
3216 fe1d -> CalcDShape(ip, dshape1dx);
3217 fe1d -> CalcDShape(ipy, dshape1dy);
3218 fe1d -> CalcDShape(ipz, dshape1dz);
3219
3220 for (int n = 0; n < dof; n++)
3221 {
3222 dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
3223 dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
3224 dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
3225 }
3226}
3227
3229{
3230 delete fe1d;
3231
3232 delete [] I;
3233 delete [] J;
3234 delete [] K;
3235}
3236
3237
3239 : NodalFiniteElement(1, Geometry::SEGMENT, 3, 4)
3240{
3241 Nodes.IntPoint(0).x = 0.0;
3242 Nodes.IntPoint(1).x = 1.0;
3243 Nodes.IntPoint(2).x = 0.5;
3244}
3245
3247 Vector &shape) const
3248{
3249 real_t x = ip.x;
3250
3251 if (x <= 0.5)
3252 {
3253 shape(0) = 1.0 - 2.0 * x;
3254 shape(1) = 0.0;
3255 shape(2) = 2.0 * x;
3256 }
3257 else
3258 {
3259 shape(0) = 0.0;
3260 shape(1) = 2.0 * x - 1.0;
3261 shape(2) = 2.0 - 2.0 * x;
3262 }
3263}
3264
3266 DenseMatrix &dshape) const
3267{
3268 real_t x = ip.x;
3269
3270 if (x <= 0.5)
3271 {
3272 dshape(0,0) = - 2.0;
3273 dshape(1,0) = 0.0;
3274 dshape(2,0) = 2.0;
3275 }
3276 else
3277 {
3278 dshape(0,0) = 0.0;
3279 dshape(1,0) = 2.0;
3280 dshape(2,0) = - 2.0;
3281 }
3282}
3283
3285 : NodalFiniteElement(2, Geometry::TRIANGLE, 6, 5)
3286{
3287 Nodes.IntPoint(0).x = 0.0;
3288 Nodes.IntPoint(0).y = 0.0;
3289 Nodes.IntPoint(1).x = 1.0;
3290 Nodes.IntPoint(1).y = 0.0;
3291 Nodes.IntPoint(2).x = 0.0;
3292 Nodes.IntPoint(2).y = 1.0;
3293 Nodes.IntPoint(3).x = 0.5;
3294 Nodes.IntPoint(3).y = 0.0;
3295 Nodes.IntPoint(4).x = 0.5;
3296 Nodes.IntPoint(4).y = 0.5;
3297 Nodes.IntPoint(5).x = 0.0;
3298 Nodes.IntPoint(5).y = 0.5;
3299}
3300
3302 Vector &shape) const
3303{
3304 int i;
3305
3306 real_t L0, L1, L2;
3307 L0 = 2.0 * ( 1. - ip.x - ip.y );
3308 L1 = 2.0 * ( ip.x );
3309 L2 = 2.0 * ( ip.y );
3310
3311 // The reference triangle is split in 4 triangles as follows:
3312 //
3313 // T0 - 0,3,5
3314 // T1 - 1,3,4
3315 // T2 - 2,4,5
3316 // T3 - 3,4,5
3317
3318 for (i = 0; i < 6; i++)
3319 {
3320 shape(i) = 0.0;
3321 }
3322
3323 if (L0 >= 1.0) // T0
3324 {
3325 shape(0) = L0 - 1.0;
3326 shape(3) = L1;
3327 shape(5) = L2;
3328 }
3329 else if (L1 >= 1.0) // T1
3330 {
3331 shape(3) = L0;
3332 shape(1) = L1 - 1.0;
3333 shape(4) = L2;
3334 }
3335 else if (L2 >= 1.0) // T2
3336 {
3337 shape(5) = L0;
3338 shape(4) = L1;
3339 shape(2) = L2 - 1.0;
3340 }
3341 else // T3
3342 {
3343 shape(3) = 1.0 - L2;
3344 shape(4) = 1.0 - L0;
3345 shape(5) = 1.0 - L1;
3346 }
3347}
3348
3350 DenseMatrix &dshape) const
3351{
3352 int i,j;
3353
3354 real_t L0, L1, L2;
3355 L0 = 2.0 * ( 1. - ip.x - ip.y );
3356 L1 = 2.0 * ( ip.x );
3357 L2 = 2.0 * ( ip.y );
3358
3359 real_t DL0[2], DL1[2], DL2[2];
3360 DL0[0] = -2.0; DL0[1] = -2.0;
3361 DL1[0] = 2.0; DL1[1] = 0.0;
3362 DL2[0] = 0.0; DL2[1] = 2.0;
3363
3364 for (i = 0; i < 6; i++)
3365 for (j = 0; j < 2; j++)
3366 {
3367 dshape(i,j) = 0.0;
3368 }
3369
3370 if (L0 >= 1.0) // T0
3371 {
3372 for (j = 0; j < 2; j++)
3373 {
3374 dshape(0,j) = DL0[j];
3375 dshape(3,j) = DL1[j];
3376 dshape(5,j) = DL2[j];
3377 }
3378 }
3379 else if (L1 >= 1.0) // T1
3380 {
3381 for (j = 0; j < 2; j++)
3382 {
3383 dshape(3,j) = DL0[j];
3384 dshape(1,j) = DL1[j];
3385 dshape(4,j) = DL2[j];
3386 }
3387 }
3388 else if (L2 >= 1.0) // T2
3389 {
3390 for (j = 0; j < 2; j++)
3391 {
3392 dshape(5,j) = DL0[j];
3393 dshape(4,j) = DL1[j];
3394 dshape(2,j) = DL2[j];
3395 }
3396 }
3397 else // T3
3398 {
3399 for (j = 0; j < 2; j++)
3400 {
3401 dshape(3,j) = - DL2[j];
3402 dshape(4,j) = - DL0[j];
3403 dshape(5,j) = - DL1[j];
3404 }
3405 }
3406}
3407
3409 : NodalFiniteElement(3, Geometry::TETRAHEDRON, 10, 4)
3410{
3411 Nodes.IntPoint(0).x = 0.0;
3412 Nodes.IntPoint(0).y = 0.0;
3413 Nodes.IntPoint(0).z = 0.0;
3414 Nodes.IntPoint(1).x = 1.0;
3415 Nodes.IntPoint(1).y = 0.0;
3416 Nodes.IntPoint(1).z = 0.0;
3417 Nodes.IntPoint(2).x = 0.0;
3418 Nodes.IntPoint(2).y = 1.0;
3419 Nodes.IntPoint(2).z = 0.0;
3420 Nodes.IntPoint(3).x = 0.0;
3421 Nodes.IntPoint(3).y = 0.0;
3422 Nodes.IntPoint(3).z = 1.0;
3423 Nodes.IntPoint(4).x = 0.5;
3424 Nodes.IntPoint(4).y = 0.0;
3425 Nodes.IntPoint(4).z = 0.0;
3426 Nodes.IntPoint(5).x = 0.0;
3427 Nodes.IntPoint(5).y = 0.5;
3428 Nodes.IntPoint(5).z = 0.0;
3429 Nodes.IntPoint(6).x = 0.0;
3430 Nodes.IntPoint(6).y = 0.0;
3431 Nodes.IntPoint(6).z = 0.5;
3432 Nodes.IntPoint(7).x = 0.5;
3433 Nodes.IntPoint(7).y = 0.5;
3434 Nodes.IntPoint(7).z = 0.0;
3435 Nodes.IntPoint(8).x = 0.5;
3436 Nodes.IntPoint(8).y = 0.0;
3437 Nodes.IntPoint(8).z = 0.5;
3438 Nodes.IntPoint(9).x = 0.0;
3439 Nodes.IntPoint(9).y = 0.5;
3440 Nodes.IntPoint(9).z = 0.5;
3441}
3442
3444 Vector &shape) const
3445{
3446 int i;
3447
3448 real_t L0, L1, L2, L3, L4, L5;
3449 L0 = 2.0 * ( 1. - ip.x - ip.y - ip.z );
3450 L1 = 2.0 * ( ip.x );
3451 L2 = 2.0 * ( ip.y );
3452 L3 = 2.0 * ( ip.z );
3453 L4 = 2.0 * ( ip.x + ip.y );
3454 L5 = 2.0 * ( ip.y + ip.z );
3455
3456 // The reference tetrahedron is split in 8 tetrahedra as follows:
3457 //
3458 // T0 - 0,4,5,6
3459 // T1 - 1,4,7,8
3460 // T2 - 2,5,7,9
3461 // T3 - 3,6,8,9
3462 // T4 - 4,5,6,8
3463 // T5 - 4,5,7,8
3464 // T6 - 5,6,8,9
3465 // T7 - 5,7,8,9
3466
3467 for (i = 0; i < 10; i++)
3468 {
3469 shape(i) = 0.0;
3470 }
3471
3472 if (L0 >= 1.0) // T0
3473 {
3474 shape(0) = L0 - 1.0;
3475 shape(4) = L1;
3476 shape(5) = L2;
3477 shape(6) = L3;
3478 }
3479 else if (L1 >= 1.0) // T1
3480 {
3481 shape(4) = L0;
3482 shape(1) = L1 - 1.0;
3483 shape(7) = L2;
3484 shape(8) = L3;
3485 }
3486 else if (L2 >= 1.0) // T2
3487 {
3488 shape(5) = L0;
3489 shape(7) = L1;
3490 shape(2) = L2 - 1.0;
3491 shape(9) = L3;
3492 }
3493 else if (L3 >= 1.0) // T3
3494 {
3495 shape(6) = L0;
3496 shape(8) = L1;
3497 shape(9) = L2;
3498 shape(3) = L3 - 1.0;
3499 }
3500 else if ((L4 <= 1.0) && (L5 <= 1.0)) // T4
3501 {
3502 shape(4) = 1.0 - L5;
3503 shape(5) = L2;
3504 shape(6) = 1.0 - L4;
3505 shape(8) = 1.0 - L0;
3506 }
3507 else if ((L4 >= 1.0) && (L5 <= 1.0)) // T5
3508 {
3509 shape(4) = 1.0 - L5;
3510 shape(5) = 1.0 - L1;
3511 shape(7) = L4 - 1.0;
3512 shape(8) = L3;
3513 }
3514 else if ((L4 <= 1.0) && (L5 >= 1.0)) // T6
3515 {
3516 shape(5) = 1.0 - L3;
3517 shape(6) = 1.0 - L4;
3518 shape(8) = L1;
3519 shape(9) = L5 - 1.0;
3520 }
3521 else if ((L4 >= 1.0) && (L5 >= 1.0)) // T7
3522 {
3523 shape(5) = L0;
3524 shape(7) = L4 - 1.0;
3525 shape(8) = 1.0 - L2;
3526 shape(9) = L5 - 1.0;
3527 }
3528}
3529
3531 DenseMatrix &dshape) const
3532{
3533 int i,j;
3534
3535 real_t L0, L1, L2, L3, L4, L5;
3536 L0 = 2.0 * ( 1. - ip.x - ip.y - ip.z );
3537 L1 = 2.0 * ( ip.x );
3538 L2 = 2.0 * ( ip.y );
3539 L3 = 2.0 * ( ip.z );
3540 L4 = 2.0 * ( ip.x + ip.y );
3541 L5 = 2.0 * ( ip.y + ip.z );
3542
3543 real_t DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
3544 DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
3545 DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
3546 DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
3547 DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
3548 DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
3549 DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
3550
3551 for (i = 0; i < 10; i++)
3552 for (j = 0; j < 3; j++)
3553 {
3554 dshape(i,j) = 0.0;
3555 }
3556
3557 if (L0 >= 1.0) // T0
3558 {
3559 for (j = 0; j < 3; j++)
3560 {
3561 dshape(0,j) = DL0[j];
3562 dshape(4,j) = DL1[j];
3563 dshape(5,j) = DL2[j];
3564 dshape(6,j) = DL3[j];
3565 }
3566 }
3567 else if (L1 >= 1.0) // T1
3568 {
3569 for (j = 0; j < 3; j++)
3570 {
3571 dshape(4,j) = DL0[j];
3572 dshape(1,j) = DL1[j];
3573 dshape(7,j) = DL2[j];
3574 dshape(8,j) = DL3[j];
3575 }
3576 }
3577 else if (L2 >= 1.0) // T2
3578 {
3579 for (j = 0; j < 3; j++)
3580 {
3581 dshape(5,j) = DL0[j];
3582 dshape(7,j) = DL1[j];
3583 dshape(2,j) = DL2[j];
3584 dshape(9,j) = DL3[j];
3585 }
3586 }
3587 else if (L3 >= 1.0) // T3
3588 {
3589 for (j = 0; j < 3; j++)
3590 {
3591 dshape(6,j) = DL0[j];
3592 dshape(8,j) = DL1[j];
3593 dshape(9,j) = DL2[j];
3594 dshape(3,j) = DL3[j];
3595 }
3596 }
3597 else if ((L4 <= 1.0) && (L5 <= 1.0)) // T4
3598 {
3599 for (j = 0; j < 3; j++)
3600 {
3601 dshape(4,j) = - DL5[j];
3602 dshape(5,j) = DL2[j];
3603 dshape(6,j) = - DL4[j];
3604 dshape(8,j) = - DL0[j];
3605 }
3606 }
3607 else if ((L4 >= 1.0) && (L5 <= 1.0)) // T5
3608 {
3609 for (j = 0; j < 3; j++)
3610 {
3611 dshape(4,j) = - DL5[j];
3612 dshape(5,j) = - DL1[j];
3613 dshape(7,j) = DL4[j];
3614 dshape(8,j) = DL3[j];
3615 }
3616 }
3617 else if ((L4 <= 1.0) && (L5 >= 1.0)) // T6
3618 {
3619 for (j = 0; j < 3; j++)
3620 {
3621 dshape(5,j) = - DL3[j];
3622 dshape(6,j) = - DL4[j];
3623 dshape(8,j) = DL1[j];
3624 dshape(9,j) = DL5[j];
3625 }
3626 }
3627 else if ((L4 >= 1.0) && (L5 >= 1.0)) // T7
3628 {
3629 for (j = 0; j < 3; j++)
3630 {
3631 dshape(5,j) = DL0[j];
3632 dshape(7,j) = DL4[j];
3633 dshape(8,j) = - DL2[j];
3634 dshape(9,j) = DL5[j];
3635 }
3636 }
3637}
3638
3639
3641 : NodalFiniteElement(2, Geometry::SQUARE, 9, 1, FunctionSpace::rQk)
3642{
3643 Nodes.IntPoint(0).x = 0.0;
3644 Nodes.IntPoint(0).y = 0.0;
3645 Nodes.IntPoint(1).x = 1.0;
3646 Nodes.IntPoint(1).y = 0.0;
3647 Nodes.IntPoint(2).x = 1.0;
3648 Nodes.IntPoint(2).y = 1.0;
3649 Nodes.IntPoint(3).x = 0.0;
3650 Nodes.IntPoint(3).y = 1.0;
3651 Nodes.IntPoint(4).x = 0.5;
3652 Nodes.IntPoint(4).y = 0.0;
3653 Nodes.IntPoint(5).x = 1.0;
3654 Nodes.IntPoint(5).y = 0.5;
3655 Nodes.IntPoint(6).x = 0.5;
3656 Nodes.IntPoint(6).y = 1.0;
3657 Nodes.IntPoint(7).x = 0.0;
3658 Nodes.IntPoint(7).y = 0.5;
3659 Nodes.IntPoint(8).x = 0.5;
3660 Nodes.IntPoint(8).y = 0.5;
3661}
3662
3664 Vector &shape) const
3665{
3666 int i;
3667 real_t x = ip.x, y = ip.y;
3668 real_t Lx, Ly;
3669 Lx = 2.0 * ( 1. - x );
3670 Ly = 2.0 * ( 1. - y );
3671
3672 // The reference square is split in 4 squares as follows:
3673 //
3674 // T0 - 0,4,7,8
3675 // T1 - 1,4,5,8
3676 // T2 - 2,5,6,8
3677 // T3 - 3,6,7,8
3678
3679 for (i = 0; i < 9; i++)
3680 {
3681 shape(i) = 0.0;
3682 }
3683
3684 if ((x <= 0.5) && (y <= 0.5)) // T0
3685 {
3686 shape(0) = (Lx - 1.0) * (Ly - 1.0);
3687 shape(4) = (2.0 - Lx) * (Ly - 1.0);
3688 shape(8) = (2.0 - Lx) * (2.0 - Ly);
3689 shape(7) = (Lx - 1.0) * (2.0 - Ly);
3690 }
3691 else if ((x >= 0.5) && (y <= 0.5)) // T1
3692 {
3693 shape(4) = Lx * (Ly - 1.0);
3694 shape(1) = (1.0 - Lx) * (Ly - 1.0);
3695 shape(5) = (1.0 - Lx) * (2.0 - Ly);
3696 shape(8) = Lx * (2.0 - Ly);
3697 }
3698 else if ((x >= 0.5) && (y >= 0.5)) // T2
3699 {
3700 shape(8) = Lx * Ly ;
3701 shape(5) = (1.0 - Lx) * Ly ;
3702 shape(2) = (1.0 - Lx) * (1.0 - Ly);
3703 shape(6) = Lx * (1.0 - Ly);
3704 }
3705 else if ((x <= 0.5) && (y >= 0.5)) // T3
3706 {
3707 shape(7) = (Lx - 1.0) * Ly ;
3708 shape(8) = (2.0 - Lx) * Ly ;
3709 shape(6) = (2.0 - Lx) * (1.0 - Ly);
3710 shape(3) = (Lx - 1.0) * (1.0 - Ly);
3711 }
3712}
3713
3715 DenseMatrix &dshape) const
3716{
3717 int i,j;
3718 real_t x = ip.x, y = ip.y;
3719 real_t Lx, Ly;
3720 Lx = 2.0 * ( 1. - x );
3721 Ly = 2.0 * ( 1. - y );
3722
3723 for (i = 0; i < 9; i++)
3724 for (j = 0; j < 2; j++)
3725 {
3726 dshape(i,j) = 0.0;
3727 }
3728
3729 if ((x <= 0.5) && (y <= 0.5)) // T0
3730 {
3731 dshape(0,0) = 2.0 * (1.0 - Ly);
3732 dshape(0,1) = 2.0 * (1.0 - Lx);
3733
3734 dshape(4,0) = 2.0 * (Ly - 1.0);
3735 dshape(4,1) = -2.0 * (2.0 - Lx);
3736
3737 dshape(8,0) = 2.0 * (2.0 - Ly);
3738 dshape(8,1) = 2.0 * (2.0 - Lx);
3739
3740 dshape(7,0) = -2.0 * (2.0 - Ly);
3741 dshape(7,0) = 2.0 * (Lx - 1.0);
3742 }
3743 else if ((x >= 0.5) && (y <= 0.5)) // T1
3744 {
3745 dshape(4,0) = -2.0 * (Ly - 1.0);
3746 dshape(4,1) = -2.0 * Lx;
3747
3748 dshape(1,0) = 2.0 * (Ly - 1.0);
3749 dshape(1,1) = -2.0 * (1.0 - Lx);
3750
3751 dshape(5,0) = 2.0 * (2.0 - Ly);
3752 dshape(5,1) = 2.0 * (1.0 - Lx);
3753
3754 dshape(8,0) = -2.0 * (2.0 - Ly);
3755 dshape(8,1) = 2.0 * Lx;
3756 }
3757 else if ((x >= 0.5) && (y >= 0.5)) // T2
3758 {
3759 dshape(8,0) = -2.0 * Ly;
3760 dshape(8,1) = -2.0 * Lx;
3761
3762 dshape(5,0) = 2.0 * Ly;
3763 dshape(5,1) = -2.0 * (1.0 - Lx);
3764
3765 dshape(2,0) = 2.0 * (1.0 - Ly);
3766 dshape(2,1) = 2.0 * (1.0 - Lx);
3767
3768 dshape(6,0) = -2.0 * (1.0 - Ly);
3769 dshape(6,1) = 2.0 * Lx;
3770 }
3771 else if ((x <= 0.5) && (y >= 0.5)) // T3
3772 {
3773 dshape(7,0) = -2.0 * Ly;
3774 dshape(7,1) = -2.0 * (Lx - 1.0);
3775
3776 dshape(8,0) = 2.0 * Ly ;
3777 dshape(8,1) = -2.0 * (2.0 - Lx);
3778
3779 dshape(6,0) = 2.0 * (1.0 - Ly);
3780 dshape(6,1) = 2.0 * (2.0 - Lx);
3781
3782 dshape(3,0) = -2.0 * (1.0 - Ly);
3783 dshape(3,1) = 2.0 * (Lx - 1.0);
3784 }
3785}
3786
3788 : NodalFiniteElement(3, Geometry::CUBE, 27, 2, FunctionSpace::rQk)
3789{
3790 real_t I[27];
3791 real_t J[27];
3792 real_t K[27];
3793 // nodes
3794 I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
3795 I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
3796 I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
3797 I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
3798 I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
3799 I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
3800 I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
3801 I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
3802 // edges
3803 I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
3804 I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
3805 I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
3806 I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
3807 I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
3808 I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
3809 I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
3810 I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
3811 I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
3812 I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
3813 I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
3814 I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
3815 // faces
3816 I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
3817 I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
3818 I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
3819 I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
3820 I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
3821 I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
3822 // element
3823 I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
3824
3825 for (int n = 0; n < 27; n++)
3826 {
3827 Nodes.IntPoint(n).x = I[n];
3828 Nodes.IntPoint(n).y = J[n];
3829 Nodes.IntPoint(n).z = K[n];
3830 }
3831}
3832
3834 Vector &shape) const
3835{
3836 int i, N[8];
3837 real_t Lx, Ly, Lz;
3838 real_t x = ip.x, y = ip.y, z = ip.z;
3839
3840 for (i = 0; i < 27; i++)
3841 {
3842 shape(i) = 0.0;
3843 }
3844
3845 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5)) // T0
3846 {
3847 Lx = 1.0 - 2.0 * x;
3848 Ly = 1.0 - 2.0 * y;
3849 Lz = 1.0 - 2.0 * z;
3850
3851 N[0] = 0;
3852 N[1] = 8;
3853 N[2] = 20;
3854 N[3] = 11;
3855 N[4] = 16;
3856 N[5] = 21;
3857 N[6] = 26;
3858 N[7] = 24;
3859 }
3860 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5)) // T1
3861 {
3862 Lx = 2.0 - 2.0 * x;
3863 Ly = 1.0 - 2.0 * y;
3864 Lz = 1.0 - 2.0 * z;
3865
3866 N[0] = 8;
3867 N[1] = 1;
3868 N[2] = 9;
3869 N[3] = 20;
3870 N[4] = 21;
3871 N[5] = 17;
3872 N[6] = 22;
3873 N[7] = 26;
3874 }
3875 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5)) // T2
3876 {
3877 Lx = 2.0 - 2.0 * x;
3878 Ly = 2.0 - 2.0 * y;
3879 Lz = 1.0 - 2.0 * z;
3880
3881 N[0] = 20;
3882 N[1] = 9;
3883 N[2] = 2;
3884 N[3] = 10;
3885 N[4] = 26;
3886 N[5] = 22;
3887 N[6] = 18;
3888 N[7] = 23;
3889 }
3890 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5)) // T3
3891 {
3892 Lx = 1.0 - 2.0 * x;
3893 Ly = 2.0 - 2.0 * y;
3894 Lz = 1.0 - 2.0 * z;
3895
3896 N[0] = 11;
3897 N[1] = 20;
3898 N[2] = 10;
3899 N[3] = 3;
3900 N[4] = 24;
3901 N[5] = 26;
3902 N[6] = 23;
3903 N[7] = 19;
3904 }
3905 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5)) // T4
3906 {
3907 Lx = 1.0 - 2.0 * x;
3908 Ly = 1.0 - 2.0 * y;
3909 Lz = 2.0 - 2.0 * z;
3910
3911 N[0] = 16;
3912 N[1] = 21;
3913 N[2] = 26;
3914 N[3] = 24;
3915 N[4] = 4;
3916 N[5] = 12;
3917 N[6] = 25;
3918 N[7] = 15;
3919 }
3920 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5)) // T5
3921 {
3922 Lx = 2.0 - 2.0 * x;
3923 Ly = 1.0 - 2.0 * y;
3924 Lz = 2.0 - 2.0 * z;
3925
3926 N[0] = 21;
3927 N[1] = 17;
3928 N[2] = 22;
3929 N[3] = 26;
3930 N[4] = 12;
3931 N[5] = 5;
3932 N[6] = 13;
3933 N[7] = 25;
3934 }
3935 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5)) // T6
3936 {
3937 Lx = 2.0 - 2.0 * x;
3938 Ly = 2.0 - 2.0 * y;
3939 Lz = 2.0 - 2.0 * z;
3940
3941 N[0] = 26;
3942 N[1] = 22;
3943 N[2] = 18;
3944 N[3] = 23;
3945 N[4] = 25;
3946 N[5] = 13;
3947 N[6] = 6;
3948 N[7] = 14;
3949 }
3950 else // T7
3951 {
3952 Lx = 1.0 - 2.0 * x;
3953 Ly = 2.0 - 2.0 * y;
3954 Lz = 2.0 - 2.0 * z;
3955
3956 N[0] = 24;
3957 N[1] = 26;
3958 N[2] = 23;
3959 N[3] = 19;
3960 N[4] = 15;
3961 N[5] = 25;
3962 N[6] = 14;
3963 N[7] = 7;
3964 }
3965
3966 shape(N[0]) = Lx * Ly * Lz;
3967 shape(N[1]) = (1 - Lx) * Ly * Lz;
3968 shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
3969 shape(N[3]) = Lx * (1 - Ly) * Lz;
3970 shape(N[4]) = Lx * Ly * (1 - Lz);
3971 shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
3972 shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
3973 shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
3974}
3975
3977 DenseMatrix &dshape) const
3978{
3979 int i, j, N[8];
3980 real_t Lx, Ly, Lz;
3981 real_t x = ip.x, y = ip.y, z = ip.z;
3982
3983 for (i = 0; i < 27; i++)
3984 for (j = 0; j < 3; j++)
3985 {
3986 dshape(i,j) = 0.0;
3987 }
3988
3989 if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5)) // T0
3990 {
3991 Lx = 1.0 - 2.0 * x;
3992 Ly = 1.0 - 2.0 * y;
3993 Lz = 1.0 - 2.0 * z;
3994
3995 N[0] = 0;
3996 N[1] = 8;
3997 N[2] = 20;
3998 N[3] = 11;
3999 N[4] = 16;
4000 N[5] = 21;
4001 N[6] = 26;
4002 N[7] = 24;
4003 }
4004 else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5)) // T1
4005 {
4006 Lx = 2.0 - 2.0 * x;
4007 Ly = 1.0 - 2.0 * y;
4008 Lz = 1.0 - 2.0 * z;
4009
4010 N[0] = 8;
4011 N[1] = 1;
4012 N[2] = 9;
4013 N[3] = 20;
4014 N[4] = 21;
4015 N[5] = 17;
4016 N[6] = 22;
4017 N[7] = 26;
4018 }
4019 else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5)) // T2
4020 {
4021 Lx = 2.0 - 2.0 * x;
4022 Ly = 2.0 - 2.0 * y;
4023 Lz = 1.0 - 2.0 * z;
4024
4025 N[0] = 20;
4026 N[1] = 9;
4027 N[2] = 2;
4028 N[3] = 10;
4029 N[4] = 26;
4030 N[5] = 22;
4031 N[6] = 18;
4032 N[7] = 23;
4033 }
4034 else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5)) // T3
4035 {
4036 Lx = 1.0 - 2.0 * x;
4037 Ly = 2.0 - 2.0 * y;
4038 Lz = 1.0 - 2.0 * z;
4039
4040 N[0] = 11;
4041 N[1] = 20;
4042 N[2] = 10;
4043 N[3] = 3;
4044 N[4] = 24;
4045 N[5] = 26;
4046 N[6] = 23;
4047 N[7] = 19;
4048 }
4049 else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5)) // T4
4050 {
4051 Lx = 1.0 - 2.0 * x;
4052 Ly = 1.0 - 2.0 * y;
4053 Lz = 2.0 - 2.0 * z;
4054
4055 N[0] = 16;
4056 N[1] = 21;
4057 N[2] = 26;
4058 N[3] = 24;
4059 N[4] = 4;
4060 N[5] = 12;
4061 N[6] = 25;
4062 N[7] = 15;
4063 }
4064 else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5)) // T5
4065 {
4066 Lx = 2.0 - 2.0 * x;
4067 Ly = 1.0 - 2.0 * y;
4068 Lz = 2.0 - 2.0 * z;
4069
4070 N[0] = 21;
4071 N[1] = 17;
4072 N[2] = 22;
4073 N[3] = 26;
4074 N[4] = 12;
4075 N[5] = 5;
4076 N[6] = 13;
4077 N[7] = 25;
4078 }
4079 else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5)) // T6
4080 {
4081 Lx = 2.0 - 2.0 * x;
4082 Ly = 2.0 - 2.0 * y;
4083 Lz = 2.0 - 2.0 * z;
4084
4085 N[0] = 26;
4086 N[1] = 22;
4087 N[2] = 18;
4088 N[3] = 23;
4089 N[4] = 25;
4090 N[5] = 13;
4091 N[6] = 6;
4092 N[7] = 14;
4093 }
4094 else // T7
4095 {
4096 Lx = 1.0 - 2.0 * x;
4097 Ly = 2.0 - 2.0 * y;
4098 Lz = 2.0 - 2.0 * z;
4099
4100 N[0] = 24;
4101 N[1] = 26;
4102 N[2] = 23;
4103 N[3] = 19;
4104 N[4] = 15;
4105 N[5] = 25;
4106 N[6] = 14;
4107 N[7] = 7;
4108 }
4109
4110 dshape(N[0],0) = -2.0 * Ly * Lz ;
4111 dshape(N[0],1) = -2.0 * Lx * Lz ;
4112 dshape(N[0],2) = -2.0 * Lx * Ly ;
4113
4114 dshape(N[1],0) = 2.0 * Ly * Lz ;
4115 dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4116 dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4117
4118 dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4119 dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4120 dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4121
4122 dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4123 dshape(N[3],1) = 2.0 * Lx * Lz ;
4124 dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4125
4126 dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4127 dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4128 dshape(N[4],2) = 2.0 * Lx * Ly ;
4129
4130 dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4131 dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4132 dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4133
4134 dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4135 dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4136 dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
4137
4138 dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
4139 dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
4140 dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
4141}
4142
4143
4145 : VectorFiniteElement(3, Geometry::CUBE, 12, 1, H_CURL,
4146 FunctionSpace::Qk)
4147{
4148 // not real nodes ...
4149 Nodes.IntPoint(0).x = 0.5;
4150 Nodes.IntPoint(0).y = 0.0;
4151 Nodes.IntPoint(0).z = 0.0;
4152
4153 Nodes.IntPoint(1).x = 1.0;
4154 Nodes.IntPoint(1).y = 0.5;
4155 Nodes.IntPoint(1).z = 0.0;
4156
4157 Nodes.IntPoint(2).x = 0.5;
4158 Nodes.IntPoint(2).y = 1.0;
4159 Nodes.IntPoint(2).z = 0.0;
4160
4161 Nodes.IntPoint(3).x = 0.0;
4162 Nodes.IntPoint(3).y = 0.5;
4163 Nodes.IntPoint(3).z = 0.0;
4164
4165 Nodes.IntPoint(4).x = 0.5;
4166 Nodes.IntPoint(4).y = 0.0;
4167 Nodes.IntPoint(4).z = 1.0;
4168
4169 Nodes.IntPoint(5).x = 1.0;
4170 Nodes.IntPoint(5).y = 0.5;
4171 Nodes.IntPoint(5).z = 1.0;
4172
4173 Nodes.IntPoint(6).x = 0.5;
4174 Nodes.IntPoint(6).y = 1.0;
4175 Nodes.IntPoint(6).z = 1.0;
4176
4177 Nodes.IntPoint(7).x = 0.0;
4178 Nodes.IntPoint(7).y = 0.5;
4179 Nodes.IntPoint(7).z = 1.0;
4180
4181 Nodes.IntPoint(8).x = 0.0;
4182 Nodes.IntPoint(8).y = 0.0;
4183 Nodes.IntPoint(8).z = 0.5;
4184
4185 Nodes.IntPoint(9).x = 1.0;
4186 Nodes.IntPoint(9).y = 0.0;
4187 Nodes.IntPoint(9).z = 0.5;
4188
4189 Nodes.IntPoint(10).x= 1.0;
4190 Nodes.IntPoint(10).y= 1.0;
4191 Nodes.IntPoint(10).z= 0.5;
4192
4193 Nodes.IntPoint(11).x= 0.0;
4194 Nodes.IntPoint(11).y= 1.0;
4195 Nodes.IntPoint(11).z= 0.5;
4196}
4197
4199 DenseMatrix &shape) const
4200{
4201 real_t x = ip.x, y = ip.y, z = ip.z;
4202
4203 shape(0,0) = (1. - y) * (1. - z);
4204 shape(0,1) = 0.;
4205 shape(0,2) = 0.;
4206
4207 shape(2,0) = y * (1. - z);
4208 shape(2,1) = 0.;
4209 shape(2,2) = 0.;
4210
4211 shape(4,0) = z * (1. - y);
4212 shape(4,1) = 0.;
4213 shape(4,2) = 0.;
4214
4215 shape(6,0) = y * z;
4216 shape(6,1) = 0.;
4217 shape(6,2) = 0.;
4218
4219 shape(1,0) = 0.;
4220 shape(1,1) = x * (1. - z);
4221 shape(1,2) = 0.;
4222
4223 shape(3,0) = 0.;
4224 shape(3,1) = (1. - x) * (1. - z);
4225 shape(3,2) = 0.;
4226
4227 shape(5,0) = 0.;
4228 shape(5,1) = x * z;
4229 shape(5,2) = 0.;
4230
4231 shape(7,0) = 0.;
4232 shape(7,1) = (1. - x) * z;
4233 shape(7,2) = 0.;
4234
4235 shape(8,0) = 0.;
4236 shape(8,1) = 0.;
4237 shape(8,2) = (1. - x) * (1. - y);
4238
4239 shape(9,0) = 0.;
4240 shape(9,1) = 0.;
4241 shape(9,2) = x * (1. - y);
4242
4243 shape(10,0) = 0.;
4244 shape(10,1) = 0.;
4245 shape(10,2) = x * y;
4246
4247 shape(11,0) = 0.;
4248 shape(11,1) = 0.;
4249 shape(11,2) = y * (1. - x);
4250
4251}
4252
4254 DenseMatrix &curl_shape)
4255const
4256{
4257 real_t x = ip.x, y = ip.y, z = ip.z;
4258
4259 curl_shape(0,0) = 0.;
4260 curl_shape(0,1) = y - 1.;
4261 curl_shape(0,2) = 1. - z;
4262
4263 curl_shape(2,0) = 0.;
4264 curl_shape(2,1) = -y;
4265 curl_shape(2,2) = z - 1.;
4266
4267 curl_shape(4,0) = 0;
4268 curl_shape(4,1) = 1. - y;
4269 curl_shape(4,2) = z;
4270
4271 curl_shape(6,0) = 0.;
4272 curl_shape(6,1) = y;
4273 curl_shape(6,2) = -z;
4274
4275 curl_shape(1,0) = x;
4276 curl_shape(1,1) = 0.;
4277 curl_shape(1,2) = 1. - z;
4278
4279 curl_shape(3,0) = 1. - x;
4280 curl_shape(3,1) = 0.;
4281 curl_shape(3,2) = z - 1.;
4282
4283 curl_shape(5,0) = -x;
4284 curl_shape(5,1) = 0.;
4285 curl_shape(5,2) = z;
4286
4287 curl_shape(7,0) = x - 1.;
4288 curl_shape(7,1) = 0.;
4289 curl_shape(7,2) = -z;
4290
4291 curl_shape(8,0) = x - 1.;
4292 curl_shape(8,1) = 1. - y;
4293 curl_shape(8,2) = 0.;
4294
4295 curl_shape(9,0) = -x;
4296 curl_shape(9,1) = y - 1.;
4297 curl_shape(9,2) = 0;
4298
4299 curl_shape(10,0) = x;
4300 curl_shape(10,1) = -y;
4301 curl_shape(10,2) = 0.;
4302
4303 curl_shape(11,0) = 1. - x;
4304 curl_shape(11,1) = y;
4305 curl_shape(11,2) = 0.;
4306}
4307
4308const real_t Nedelec1HexFiniteElement::tk[12][3] =
4309{
4310 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4311 {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4312 {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
4313};
4314
4316 ElementTransformation &Trans, DenseMatrix &I) const
4317{
4318 int k, j;
4321#endif
4322
4323#ifdef MFEM_DEBUG
4324 for (k = 0; k < dof; k++)
4325 {
4327 for (j = 0; j < dof; j++)
4328 {
4329 real_t d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
4330 vshape(j,2)*tk[k][2] );
4331 if (j == k) { d -= 1.0; }
4332 if (fabs(d) > 1.0e-12)
4333 {
4334 mfem::err << "Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
4335 " k = " << k << ", j = " << j << ", d = " << d << endl;
4336 mfem_error();
4337 }
4338 }
4339 }
4340#endif
4341
4343 ip.x = ip.y = ip.z = 0.0;
4344 Trans.SetIntPoint (&ip);
4345 // Trans must be linear (more to have embedding?)
4346 const DenseMatrix &J = Trans.Jacobian();
4347 real_t vk[3];
4348 Vector xk (vk, 3);
4349
4350 for (k = 0; k < dof; k++)
4351 {
4352 Trans.Transform (Nodes.IntPoint (k), xk);
4353 ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
4354 CalcVShape (ip, vshape);
4355 // vk = J tk
4356 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4357 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4358 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4359 for (j = 0; j < dof; j++)
4360 if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
4361 vshape(j,2)*vk[2])) < 1.0e-12)
4362 {
4363 I(k,j) = 0.0;
4364 }
4365 }
4366}
4367
4370 Vector &dofs) const
4371{
4372 real_t vk[3];
4373 Vector xk (vk, 3);
4374
4375 for (int k = 0; k < dof; k++)
4376 {
4377 Trans.SetIntPoint (&Nodes.IntPoint (k));
4378 const DenseMatrix &J = Trans.Jacobian();
4379
4380 vc.Eval (xk, Trans, Nodes.IntPoint (k));
4381 // xk^t J tk
4382 dofs(k) =
4383 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4384 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4385 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4386 }
4387}
4388
4390 ElementTransformation &Trans,
4392{
4393 DenseMatrix dshape(fe.GetDof(), 3);
4395
4397 for (int k = 0; k < dof; k++)
4398 {
4399 fe.CalcDShape(Nodes.IntPoint(k), dshape);
4401 for (int j = 0; j < grad_k.Size(); j++)
4402 {
4404 }
4405 }
4406}
4407
4408
4410 : VectorFiniteElement(3, Geometry::TETRAHEDRON, 6, 1, H_CURL)
4411{
4412 // not real nodes ...
4413 Nodes.IntPoint(0).x = 0.5;
4414 Nodes.IntPoint(0).y = 0.0;
4415 Nodes.IntPoint(0).z = 0.0;
4416
4417 Nodes.IntPoint(1).x = 0.0;
4418 Nodes.IntPoint(1).y = 0.5;
4419 Nodes.IntPoint(1).z = 0.0;
4420
4421 Nodes.IntPoint(2).x = 0.0;
4422 Nodes.IntPoint(2).y = 0.0;
4423 Nodes.IntPoint(2).z = 0.5;
4424
4425 Nodes.IntPoint(3).x = 0.5;
4426 Nodes.IntPoint(3).y = 0.5;
4427 Nodes.IntPoint(3).z = 0.0;
4428
4429 Nodes.IntPoint(4).x = 0.5;
4430 Nodes.IntPoint(4).y = 0.0;
4431 Nodes.IntPoint(4).z = 0.5;
4432
4433 Nodes.IntPoint(5).x = 0.0;
4434 Nodes.IntPoint(5).y = 0.5;
4435 Nodes.IntPoint(5).z = 0.5;
4436}
4437
4439 DenseMatrix &shape) const
4440{
4441 real_t x = ip.x, y = ip.y, z = ip.z;
4442
4443 shape(0,0) = 1. - y - z;
4444 shape(0,1) = x;
4445 shape(0,2) = x;
4446
4447 shape(1,0) = y;
4448 shape(1,1) = 1. - x - z;
4449 shape(1,2) = y;
4450
4451 shape(2,0) = z;
4452 shape(2,1) = z;
4453 shape(2,2) = 1. - x - y;
4454
4455 shape(3,0) = -y;
4456 shape(3,1) = x;
4457 shape(3,2) = 0.;
4458
4459 shape(4,0) = -z;
4460 shape(4,1) = 0.;
4461 shape(4,2) = x;
4462
4463 shape(5,0) = 0.;
4464 shape(5,1) = -z;
4465 shape(5,2) = y;
4466}
4467
4469 DenseMatrix &curl_shape)
4470const
4471{
4472 curl_shape(0,0) = 0.;
4473 curl_shape(0,1) = -2.;
4474 curl_shape(0,2) = 2.;
4475
4476 curl_shape(1,0) = 2.;
4477 curl_shape(1,1) = 0.;
4478 curl_shape(1,2) = -2.;
4479
4480 curl_shape(2,0) = -2.;
4481 curl_shape(2,1) = 2.;
4482 curl_shape(2,2) = 0.;
4483
4484 curl_shape(3,0) = 0.;
4485 curl_shape(3,1) = 0.;
4486 curl_shape(3,2) = 2.;
4487
4488 curl_shape(4,0) = 0.;
4489 curl_shape(4,1) = -2.;
4490 curl_shape(4,2) = 0.;
4491
4492 curl_shape(5,0) = 2.;
4493 curl_shape(5,1) = 0.;
4494 curl_shape(5,2) = 0.;
4495}
4496
4497const real_t Nedelec1TetFiniteElement::tk[6][3] =
4498{{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
4499
4501 ElementTransformation &Trans, DenseMatrix &I) const
4502{
4503 int k, j;
4506#endif
4507
4508#ifdef MFEM_DEBUG
4509 for (k = 0; k < dof; k++)
4510 {
4512 for (j = 0; j < dof; j++)
4513 {
4514 real_t d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
4515 vshape(j,2)*tk[k][2] );
4516 if (j == k) { d -= 1.0; }
4517 if (fabs(d) > 1.0e-12)
4518 {
4519 mfem::err << "Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
4520 " k = " << k << ", j = " << j << ", d = " << d << endl;
4521 mfem_error();
4522 }
4523 }
4524 }
4525#endif
4526
4528 ip.x = ip.y = ip.z = 0.0;
4529 Trans.SetIntPoint (&ip);
4530 // Trans must be linear
4531 const DenseMatrix &J = Trans.Jacobian();
4532 real_t vk[3];
4533 Vector xk (vk, 3);
4534
4535 for (k = 0; k < dof; k++)
4536 {
4537 Trans.Transform (Nodes.IntPoint (k), xk);
4538 ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
4539 CalcVShape (ip, vshape);
4540 // vk = J tk
4541 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4542 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4543 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4544 for (j = 0; j < dof; j++)
4545 if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
4546 vshape(j,2)*vk[2])) < 1.0e-12)
4547 {
4548 I(k,j) = 0.0;
4549 }
4550 }
4551}
4552
4555 Vector &dofs) const
4556{
4557 real_t vk[3];
4558 Vector xk (vk, 3);
4559
4560 for (int k = 0; k < dof; k++)
4561 {
4562 Trans.SetIntPoint (&Nodes.IntPoint (k));
4563 const DenseMatrix &J = Trans.Jacobian();
4564
4565 vc.Eval (xk, Trans, Nodes.IntPoint (k));
4566 // xk^t J tk
4567 dofs(k) =
4568 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4569 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4570 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4571 }
4572}
4573
4575 ElementTransformation &Trans,
4577{
4578 DenseMatrix dshape(fe.GetDof(), 3);
4580
4582 for (int k = 0; k < dof; k++)
4583 {
4584 fe.CalcDShape(Nodes.IntPoint(k), dshape);
4586 for (int j = 0; j < grad_k.Size(); j++)
4587 {
4589 }
4590 }
4591}
4592
4593
4595 : VectorFiniteElement(3, Geometry::PRISM, 9, 1, H_CURL, FunctionSpace::Qk)
4596{
4597 // not real nodes ...
4598 Nodes.IntPoint(0).x = 0.5;
4599 Nodes.IntPoint(0).y = 0.0;
4600 Nodes.IntPoint(0).z = 0.0;
4601
4602 Nodes.IntPoint(1).x = 0.5;
4603 Nodes.IntPoint(1).y = 0.5;
4604 Nodes.IntPoint(1).z = 0.0;
4605
4606 Nodes.IntPoint(2).x = 0.0;
4607 Nodes.IntPoint(2).y = 0.5;
4608 Nodes.IntPoint(2).z = 0.0;
4609
4610 Nodes.IntPoint(3).x = 0.5;
4611 Nodes.IntPoint(3).y = 0.0;
4612 Nodes.IntPoint(3).z = 1.0;
4613
4614 Nodes.IntPoint(4).x = 0.5;
4615 Nodes.IntPoint(4).y = 0.5;
4616 Nodes.IntPoint(4).z = 1.0;
4617
4618 Nodes.IntPoint(5).x = 0.0;
4619 Nodes.IntPoint(5).y = 0.5;
4620 Nodes.IntPoint(5).z = 1.0;
4621
4622 Nodes.IntPoint(6).x = 0.0;
4623 Nodes.IntPoint(6).y = 0.0;
4624 Nodes.IntPoint(6).z = 0.5;
4625
4626 Nodes.IntPoint(7).x = 1.0;
4627 Nodes.IntPoint(7).y = 0.0;
4628 Nodes.IntPoint(7).z = 0.5;
4629
4630 Nodes.IntPoint(8).x = 0.0;
4631 Nodes.IntPoint(8).y = 1.0;
4632 Nodes.IntPoint(8).z = 0.5;
4633}
4634
4636 DenseMatrix &shape) const
4637{
4638 real_t x = ip.x, y = ip.y, z = ip.z;
4639
4640 shape(0,0) = (1. - y) * (1. - z);
4641 shape(0,1) = x * (1. - z);
4642 shape(0,2) = 0.;
4643
4644 shape(1,0) = - y * (1. - z);
4645 shape(1,1) = x * (1. - z);
4646 shape(1,2) = 0.;
4647
4648 shape(2,0) = - y * (1. - z);
4649 shape(2,1) = - (1. - x) * (1. - z);
4650 shape(2,2) = 0.;
4651
4652 shape(3,0) = (1. - y) * z;
4653 shape(3,1) = x * z;
4654 shape(3,2) = 0.;
4655
4656 shape(4,0) = - y * z;
4657 shape(4,1) = x * z;
4658 shape(4,2) = 0.;
4659
4660 shape(5,0) = - y * z;
4661 shape(5,1) = - (1. - x) * z;
4662 shape(5,2) = 0.;
4663
4664 shape(6,0) = 0.;
4665 shape(6,1) = 0.;
4666 shape(6,2) = 1. - x - y;
4667
4668 shape(7,0) = 0.;
4669 shape(7,1) = 0.;
4670 shape(7,2) = x;
4671
4672 shape(8,0) = 0.;
4673 shape(8,1) = 0.;
4674 shape(8,2) = y;
4675}
4676
4678 DenseMatrix &curl_shape)
4679const
4680{
4681 real_t x = ip.x, y = ip.y, z2 = 2. * ip.z;
4682
4683 curl_shape(0,0) = x;
4684 curl_shape(0,1) = - 1. + y;
4685 curl_shape(0,2) = 2. - z2;
4686
4687 curl_shape(1,0) = x;
4688 curl_shape(1,1) = y;
4689 curl_shape(1,2) = 2. - z2;
4690
4691 curl_shape(2,0) = - 1. + x;
4692 curl_shape(2,1) = y;
4693 curl_shape(2,2) = 2. - z2;
4694
4695 curl_shape(3,0) = - x;
4696 curl_shape(3,1) = 1. - y;
4697 curl_shape(3,2) = z2;
4698
4699 curl_shape(4,0) = - x;
4700 curl_shape(4,1) = - y;
4701 curl_shape(4,2) = z2;
4702
4703 curl_shape(5,0) = 1. - x;
4704 curl_shape(5,1) = - y;
4705 curl_shape(5,2) = z2;
4706
4707 curl_shape(6,0) = - 1.;
4708 curl_shape(6,1) = 1.;
4709 curl_shape(6,2) = 0.;
4710
4711 curl_shape(7,0) = 0.;
4712 curl_shape(7,1) = - 1.;
4713 curl_shape(7,2) = 0.;
4714
4715 curl_shape(8,0) = 1.;
4716 curl_shape(8,1) = 0.;
4717 curl_shape(8,2) = 0.;
4718}
4719
4720const real_t Nedelec1WdgFiniteElement::tk[9][3] =
4721{
4722 {1,0,0}, {-1,1,0}, {0,-1,0}, {1,0,0}, {-1,1,0}, {0,-1,0},
4723 {0,0,1}, {0,0,1}, {0,0,1}
4724};
4725
4727 ElementTransformation &Trans, DenseMatrix &I) const
4728{
4729 int k, j;
4732#endif
4733
4734#ifdef MFEM_DEBUG
4735 for (k = 0; k < dof; k++)
4736 {
4738 for (j = 0; j < dof; j++)
4739 {
4740 real_t d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
4741 vshape(j,2)*tk[k][2] );
4742 if (j == k) { d -= 1.0; }
4743 if (fabs(d) > 1.0e-12)
4744 {
4745 mfem::err << "Nedelec1WdgFiniteElement::GetLocalInterpolation (...)\n"
4746 " k = " << k << ", j = " << j << ", d = " << d << endl;
4747 mfem_error();
4748 }
4749 }
4750 }
4751#endif
4752
4754 ip.x = ip.y = ip.z = 0.0;
4755 Trans.SetIntPoint (&ip);
4756 // Trans must be linear
4757 const DenseMatrix &J = Trans.Jacobian();
4758 real_t vk[3];
4759 Vector xk (vk, 3);
4760
4761 for (k = 0; k < dof; k++)
4762 {
4763 Trans.Transform (Nodes.IntPoint (k), xk);
4764 ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
4765 CalcVShape (ip, vshape);
4766 // vk = J tk
4767 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4768 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4769 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4770 for (j = 0; j < dof; j++)
4771 if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
4772 vshape(j,2)*vk[2])) < 1.0e-12)
4773 {
4774 I(k,j) = 0.0;
4775 }
4776 }
4777}
4778
4781 Vector &dofs) const
4782{
4783 real_t vk[3];
4784 Vector xk (vk, 3);
4785
4786 for (int k = 0; k < dof; k++)
4787 {
4788 Trans.SetIntPoint (&Nodes.IntPoint (k));
4789 const DenseMatrix &J = Trans.Jacobian();
4790
4791 vc.Eval (xk, Trans, Nodes.IntPoint (k));
4792 // xk^t J tk
4793 dofs(k) =
4794 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4795 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4796 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4797 }
4798}
4799
4801 ElementTransformation &Trans,
4803{
4804 DenseMatrix dshape(fe.GetDof(), 3);
4806
4808 for (int k = 0; k < dof; k++)
4809 {
4810 fe.CalcDShape(Nodes.IntPoint(k), dshape);
4812 for (int j = 0; j < grad_k.Size(); j++)
4813 {
4815 }
4816 }
4817}
4818
4819
4821 : VectorFiniteElement(3, Geometry::PYRAMID, 8, 1, H_CURL)
4822{
4823 // not real nodes ...
4824 Nodes.IntPoint(0).x = 0.5;
4825 Nodes.IntPoint(0).y = 0.0;
4826 Nodes.IntPoint(0).z = 0.0;
4827
4828 Nodes.IntPoint(1).x = 1.0;
4829 Nodes.IntPoint(1).y = 0.5;
4830 Nodes.IntPoint(1).z = 0.0;
4831
4832 Nodes.IntPoint(2).x = 0.5;
4833 Nodes.IntPoint(2).y = 1.0;
4834 Nodes.IntPoint(2).z = 0.0;
4835
4836 Nodes.IntPoint(3).x = 0.0;
4837 Nodes.IntPoint(3).y = 0.5;
4838 Nodes.IntPoint(3).z = 0.0;
4839
4840 Nodes.IntPoint(4).x = 0.0;
4841 Nodes.IntPoint(4).y = 0.0;
4842 Nodes.IntPoint(4).z = 0.5;
4843
4844 Nodes.IntPoint(5).x = 0.5;
4845 Nodes.IntPoint(5).y = 0.0;
4846 Nodes.IntPoint(5).z = 0.5;
4847
4848 Nodes.IntPoint(6).x = 0.5;
4849 Nodes.IntPoint(6).y = 0.5;
4850 Nodes.IntPoint(6).z = 0.5;
4851
4852 Nodes.IntPoint(7).x = 0.0;
4853 Nodes.IntPoint(7).y = 0.5;
4854 Nodes.IntPoint(7).z = 0.5;
4855}
4856
4858 DenseMatrix &shape) const
4859{
4860 real_t x = ip.x, y = ip.y, z = ip.z, z2 = 2. * ip.z;
4861 real_t ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4862
4863 real_t tol = 1e-6;
4864
4865 if (oz <= tol)
4866 {
4867 // We must return the limit of the basis functions as z->1. In order to
4868 // remain inside the pyramid in this limit the x and y coordinates must
4869 // be approaching 0. The resulting limiting basis function values are:
4870 shape(0,0) = 0.;
4871 shape(0,1) = 0.;
4872 shape(0,2) = 0.;
4873
4874 shape(1,0) = 0.;
4875 shape(1,1) = 0.;
4876 shape(1,2) = 0.;
4877
4878 shape(2,0) = 0.;
4879 shape(2,1) = 0.;
4880 shape(2,2) = 0.;
4881
4882 shape(3,0) = 0.;
4883 shape(3,1) = 0.;
4884 shape(3,2) = 0.;
4885
4886 shape(4,0) = 1.;
4887 shape(4,1) = 1.;
4888 shape(4,2) = 1.;
4889
4890 shape(5,0) = - 1.;
4891 shape(5,1) = 0.;
4892 shape(5,2) = 0.;
4893
4894 shape(6,0) = 0.;
4895 shape(6,1) = 0.;
4896 shape(6,2) = 0.;
4897
4898 shape(7,0) = 0.;
4899 shape(7,1) = - 1.;
4900 shape(7,2) = 0.;
4901
4902 return;
4903 }
4904
4905 real_t ozi = 1.0 / oz;
4906
4907 shape(0,0) = oy;
4908 shape(0,1) = 0.;
4909 shape(0,2) = x * oy * ozi;
4910
4911 shape(1,0) = 0.;
4912 shape(1,1) = x;
4913 shape(1,2) = x * y * ozi;
4914
4915 shape(2,0) = y;
4916 shape(2,1) = 0.;
4917 shape(2,2) = x * y * ozi;
4918
4919 shape(3,0) = 0.;
4920 shape(3,1) = ox;
4921 shape(3,2) = ox * y * ozi;
4922
4923 shape(4,0) = oy * z * ozi;
4924 shape(4,1) = ox * z * ozi;
4925 shape(4,2) = 1. - x - y + x * y * (1. - z2) * ozi * ozi;
4926
4927 shape(5,0) = - oy * z * ozi;
4928 shape(5,1) = x * z * ozi;
4929 shape(5,2) = x * (1. - y * (1. - z2) * ozi * ozi);
4930
4931 shape(6,0) = - y * z * ozi;
4932 shape(6,1) = - x * z * ozi;
4933 shape(6,2) = x * y * (1. - z2) * ozi * ozi;
4934
4935 shape(7,0) = y * z * ozi;
4936 shape(7,1) = - ox * z * ozi;
4937 shape(7,2) = y * (1. - x * (1. - z2) * ozi * ozi);
4938}
4939
4941 DenseMatrix &curl_shape)
4942const
4943{
4944 real_t x = ip.x, y = ip.y, z = ip.z, z2 = 2. * z;
4945 real_t ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4946
4947 real_t tol = 1e-6;
4948
4949 if (oz <= tol)
4950 {
4951 // We must return the limit of the basis function derivatives as z->1.
4952 // In order to remain inside the pyramid in this limit the x and y
4953 // coordinates must be approaching 0. The resulting limiting basis
4954 // function values are:
4955 curl_shape(0,0) = 0.;
4956 curl_shape(0,1) = - 2.;
4957 curl_shape(0,2) = 1.;
4958
4959 curl_shape(1,0) = 0.;
4960 curl_shape(1,1) = 0.;
4961 curl_shape(1,2) = 1.;
4962
4963 curl_shape(2,0) = 0.;
4964 curl_shape(2,1) = 0.;
4965 curl_shape(2,2) = - 1.;
4966
4967 curl_shape(3,0) = 2.;
4968 curl_shape(3,1) = 0.;
4969 curl_shape(3,2) = - 1.;
4970
4971 curl_shape(4,0) = - 2.;
4972 curl_shape(4,1) = 2.;
4973 curl_shape(4,2) = 0.;
4974
4975 curl_shape(5,0) = 0.;
4976 curl_shape(5,1) = - 2.;
4977 curl_shape(5,2) = 0.;
4978
4979 curl_shape(6,0) = 0.;
4980 curl_shape(6,1) = 0.;
4981 curl_shape(6,2) = 0.;
4982
4983 curl_shape(7,0) = 2.;
4984 curl_shape(7,1) = 0.;
4985 curl_shape(7,2) = 0.;
4986
4987 return;
4988 }
4989
4990 real_t ozi = 1. / oz;
4991
4992 curl_shape(0,0) = - x * ozi;
4993 curl_shape(0,1) = - 2. + y * ozi;
4994 curl_shape(0,2) = 1.;
4995
4996 curl_shape(1,0) = x * ozi;
4997 curl_shape(1,1) = - y * ozi;
4998 curl_shape(1,2) = 1.;
4999
5000 curl_shape(2,0) = x * ozi;
5001 curl_shape(2,1) = - y * ozi;
5002 curl_shape(2,2) = - 1.;
5003
5004 curl_shape(3,0) = (2. - x - z2) * ozi;
5005 curl_shape(3,1) = y * ozi;
5006 curl_shape(3,2) = - 1.;
5007
5008 curl_shape(4,0) = - 2. * ox * ozi;
5009 curl_shape(4,1) = 2. * oy * ozi;
5010 curl_shape(4,2) = 0.;
5011
5012 curl_shape(5,0) = - 2. * x * ozi;
5013 curl_shape(5,1) = - 2. * oy * ozi;
5014 curl_shape(5,2) = 0.;
5015
5016 curl_shape(6,0) = 2. * x * ozi;
5017 curl_shape(6,1) = - 2. * y * ozi;
5018 curl_shape(6,2) = 0.;
5019
5020 curl_shape(7,0) = 2. * ox * ozi;
5021 curl_shape(7,1) = 2. * y * ozi;
5022 curl_shape(7,2) = 0.;
5023}
5024
5025const real_t Nedelec1PyrFiniteElement::tk[8][3] =
5026{{1,0,0}, {0,1,0}, {1,0,0}, {0,1,0}, {0,0,1}, {-1,0,1}, {-1,-1,1}, {0,-1,1}};
5027
5029 ElementTransformation &Trans, DenseMatrix &I) const
5030{
5031 int k, j;
5034#endif
5035
5036#ifdef MFEM_DEBUG
5037 for (k = 0; k < dof; k++)
5038 {
5040 for (j = 0; j < dof; j++)
5041 {
5042 real_t d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
5043 vshape(j,2)*tk[k][2] );
5044 if (j == k) { d -= 1.0; }
5045 if (fabs(d) > 1.0e-12)
5046 {
5047 mfem::err << "Nedelec1PyrFiniteElement::GetLocalInterpolation (...)\n"
5048 " k = " << k << ", j = " << j << ", d = " << d << endl;
5049 mfem_error();
5050 }
5051 }
5052 }
5053#endif
5054
5056 ip.x = ip.y = ip.z = 0.0;
5057 Trans.SetIntPoint (&ip);
5058 // Trans must be linear
5059 const DenseMatrix &J = Trans.Jacobian();
5060 real_t vk[3];
5061 Vector xk (vk, 3);
5062
5063 for (k = 0; k < dof; k++)
5064 {
5065 Trans.Transform (Nodes.IntPoint (k), xk);
5066 ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
5067 CalcVShape (ip, vshape);
5068 // vk = J tk
5069 vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
5070 vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
5071 vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
5072 for (j = 0; j < dof; j++)
5073 if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
5074 vshape(j,2)*vk[2])) < 1.0e-12)
5075 {
5076 I(k,j) = 0.0;
5077 }
5078 }
5079}
5080
5083 Vector &dofs) const
5084{
5085 real_t vk[3];
5086 Vector xk (vk, 3);
5087
5088 for (int k = 0; k < dof; k++)
5089 {
5090 Trans.SetIntPoint (&Nodes.IntPoint (k));
5091 const DenseMatrix &J = Trans.Jacobian();
5092
5093 vc.Eval (xk, Trans, Nodes.IntPoint (k));
5094 // xk^t J tk
5095 dofs(k) =
5096 vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
5097 vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
5098 vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
5099 }
5100}
5101
5103 ElementTransformation &Trans,
5105{
5106 DenseMatrix dshape(fe.GetDof(), 3);
5108
5110 for (int k = 0; k < dof; k++)
5111 {
5112 fe.CalcDShape(Nodes.IntPoint(k), dshape);
5114 for (int j = 0; j < grad_k.Size(); j++)
5115 {
5117 }
5118 }
5119}
5120
5121
5123 : VectorFiniteElement(3, Geometry::CUBE, 6, 1, H_DIV, FunctionSpace::Qk)
5124{
5125 // not real nodes ...
5126 // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
5127 Nodes.IntPoint(0).x = 0.5;
5128 Nodes.IntPoint(0).y = 0.5;
5129 Nodes.IntPoint(0).z = 0.0;
5130
5131 Nodes.IntPoint(1).x = 0.5;
5132 Nodes.IntPoint(1).y = 0.0;
5133 Nodes.IntPoint(1).z = 0.5;
5134
5135 Nodes.IntPoint(2).x = 1.0;
5136 Nodes.IntPoint(2).y = 0.5;
5137 Nodes.IntPoint(2).z = 0.5;
5138
5139 Nodes.IntPoint(3).x = 0.5;
5140 Nodes.IntPoint(3).y = 1.0;
5141 Nodes.IntPoint(3).