MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fe_fixed_order.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
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 
17 namespace mfem
18 {
19 
20 using namespace std;
21 
22 PointFiniteElement::PointFiniteElement()
23  : NodalFiniteElement(0, Geometry::POINT, 1, 0)
24 {
25  lex_ordering.SetSize(1);
26  lex_ordering[0] = 0;
27  Nodes.IntPoint(0).x = 0.0;
28 }
29 
30 void PointFiniteElement::CalcShape(const IntegrationPoint &ip,
31  Vector &shape) const
32 {
33  shape(0) = 1.;
34 }
35 
36 void PointFiniteElement::CalcDShape(const IntegrationPoint &ip,
37  DenseMatrix &dshape) const
38 {
39  // dshape is (1 x 0) - nothing to compute
40 }
41 
42 Linear1DFiniteElement::Linear1DFiniteElement()
43  : NodalFiniteElement(1, Geometry::SEGMENT, 2, 1)
44 {
45  Nodes.IntPoint(0).x = 0.0;
46  Nodes.IntPoint(1).x = 1.0;
47 }
48 
49 void Linear1DFiniteElement::CalcShape(const IntegrationPoint &ip,
50  Vector &shape) const
51 {
52  shape(0) = 1. - ip.x;
53  shape(1) = ip.x;
54 }
55 
56 void Linear1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
57  DenseMatrix &dshape) const
58 {
59  dshape(0,0) = -1.;
60  dshape(1,0) = 1.;
61 }
62 
63 Linear2DFiniteElement::Linear2DFiniteElement()
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 
74 void Linear2DFiniteElement::CalcShape(const IntegrationPoint &ip,
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 
82 void Linear2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
91 BiLinear2DFiniteElement::BiLinear2DFiniteElement()
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 
104 void BiLinear2DFiniteElement::CalcShape(const IntegrationPoint &ip,
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 
113 void BiLinear2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
122 void BiLinear2DFiniteElement::CalcHessian(
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 
132 GaussLinear2DFiniteElement::GaussLinear2DFiniteElement()
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 
143 void GaussLinear2DFiniteElement::CalcShape(const IntegrationPoint &ip,
144  Vector &shape) const
145 {
146  const double 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 
153 void GaussLinear2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
161 void GaussLinear2DFiniteElement::ProjectDelta(int vertex, Vector &dofs) const
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)
170 const double GaussBiLinear2DFiniteElement::p[] =
171 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
172 
173 GaussBiLinear2DFiniteElement::GaussBiLinear2DFiniteElement()
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 
186 void GaussBiLinear2DFiniteElement::CalcShape(const IntegrationPoint &ip,
187  Vector &shape) const
188 {
189  const double 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 
197 void GaussBiLinear2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
198  DenseMatrix &dshape) const
199 {
200  const double 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 
208 void GaussBiLinear2DFiniteElement::ProjectDelta(int vertex, Vector &dofs) const
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 
221 P1OnQuadFiniteElement::P1OnQuadFiniteElement()
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 
232 void P1OnQuadFiniteElement::CalcShape(const IntegrationPoint &ip,
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 
240 void P1OnQuadFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
249 Quad1DFiniteElement::Quad1DFiniteElement()
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 
257 void Quad1DFiniteElement::CalcShape(const IntegrationPoint &ip,
258  Vector &shape) const
259 {
260  double x = ip.x;
261  double 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 
268 void Quad1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
269  DenseMatrix &dshape) const
270 {
271  double 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 
279 Quad2DFiniteElement::Quad2DFiniteElement()
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 
296 void Quad2DFiniteElement::CalcShape(const IntegrationPoint &ip,
297  Vector &shape) const
298 {
299  double x = ip.x, y = ip.y;
300  double 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 
310 void Quad2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
311  DenseMatrix &dshape) const
312 {
313  double 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 
334 void Quad2DFiniteElement::CalcHessian (const IntegrationPoint &ip,
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 
362 void 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 
379 const double GaussQuad2DFiniteElement::p[] =
380 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
381 
382 GaussQuad2DFiniteElement::GaussQuad2DFiniteElement()
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 double 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 
412 void GaussQuad2DFiniteElement::CalcShape(const IntegrationPoint &ip,
413  Vector &shape) const
414 {
415  const double 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 
426 void GaussQuad2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
427  DenseMatrix &dshape) const
428 {
429  const double 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 
441 BiQuad2DFiniteElement::BiQuad2DFiniteElement()
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 
464 void BiQuad2DFiniteElement::CalcShape(const IntegrationPoint &ip,
465  Vector &shape) const
466 {
467  double x = ip.x, y = ip.y;
468  double 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 
488 void BiQuad2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
489  DenseMatrix &dshape) const
490 {
491  double x = ip.x, y = ip.y;
492  double l1x, l2x, l3x, l1y, l2y, l3y;
493  double 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 
537 void 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 
556 GaussBiQuad2DFiniteElement::GaussBiQuad2DFiniteElement()
557  : NodalFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
558 {
559  const double 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 
581 void GaussBiQuad2DFiniteElement::CalcShape(const IntegrationPoint &ip,
582  Vector &shape) const
583 {
584  const double a = sqrt(5./3.);
585  const double p1 = 0.5*(1.-sqrt(3./5.));
586 
587  double x = a*(ip.x-p1), y = a*(ip.y-p1);
588  double 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 
608 void GaussBiQuad2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
609  DenseMatrix &dshape) const
610 {
611  const double a = sqrt(5./3.);
612  const double p1 = 0.5*(1.-sqrt(3./5.));
613 
614  double x = a*(ip.x-p1), y = a*(ip.y-p1);
615  double l1x, l2x, l3x, l1y, l2y, l3y;
616  double 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 
660 BiCubic2DFiniteElement::BiCubic2DFiniteElement()
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 
697 void BiCubic2DFiniteElement::CalcShape(
698  const IntegrationPoint &ip, Vector &shape) const
699 {
700  double x = ip.x, y = ip.y;
701 
702  double w1x, w2x, w3x, w1y, w2y, w3y;
703  double 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 
736 void BiCubic2DFiniteElement::CalcDShape(
737  const IntegrationPoint &ip, DenseMatrix &dshape) const
738 {
739  double x = ip.x, y = ip.y;
740 
741  double w1x, w2x, w3x, w1y, w2y, w3y;
742  double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
743  double 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 
786 void BiCubic2DFiniteElement::CalcHessian(
787  const IntegrationPoint &ip, DenseMatrix &h) const
788 {
789  double x = ip.x, y = ip.y;
790 
791  double w1x, w2x, w3x, w1y, w2y, w3y;
792  double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
793  double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
794  double 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 
848 Cubic1DFiniteElement::Cubic1DFiniteElement()
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 
857 void Cubic1DFiniteElement::CalcShape(const IntegrationPoint &ip,
858  Vector &shape) const
859 {
860  double x = ip.x;
861  double 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 
872 void Cubic1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
873  DenseMatrix &dshape) const
874 {
875  double 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 
884 Cubic2DFiniteElement::Cubic2DFiniteElement()
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 
909 void Cubic2DFiniteElement::CalcShape(const IntegrationPoint &ip,
910  Vector &shape) const
911 {
912  double x = ip.x, y = ip.y;
913  double 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 
929 void Cubic2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
930  DenseMatrix &dshape) const
931 {
932  double 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 
957 void Cubic2DFiniteElement::CalcHessian (const IntegrationPoint &ip,
958  DenseMatrix &h) const
959 {
960  double 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 
1004 Cubic3DFiniteElement::Cubic3DFiniteElement()
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 
1069 void Cubic3DFiniteElement::CalcShape(const IntegrationPoint &ip,
1070  Vector &shape) const
1071 {
1072  double 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 
1097 void Cubic3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
1098  DenseMatrix &dshape) const
1099 {
1100  double 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 
1171 P0TriangleFiniteElement::P0TriangleFiniteElement()
1172  : NodalFiniteElement(2, Geometry::TRIANGLE, 1, 0)
1173 {
1174  Nodes.IntPoint(0).x = 0.333333333333333333;
1175  Nodes.IntPoint(0).y = 0.333333333333333333;
1176 }
1177 
1178 void P0TriangleFiniteElement::CalcShape(const IntegrationPoint &ip,
1179  Vector &shape) const
1180 {
1181  shape(0) = 1.0;
1182 }
1183 
1184 void P0TriangleFiniteElement::CalcDShape(const IntegrationPoint &ip,
1185  DenseMatrix &dshape) const
1186 {
1187  dshape(0,0) = 0.0;
1188  dshape(0,1) = 0.0;
1189 }
1190 
1191 
1192 P0QuadFiniteElement::P0QuadFiniteElement()
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 
1199 void P0QuadFiniteElement::CalcShape(const IntegrationPoint &ip,
1200  Vector &shape) const
1201 {
1202  shape(0) = 1.0;
1203 }
1204 
1205 void P0QuadFiniteElement::CalcDShape(const IntegrationPoint &ip,
1206  DenseMatrix &dshape) const
1207 {
1208  dshape(0,0) = 0.0;
1209  dshape(0,1) = 0.0;
1210 }
1211 
1212 
1213 Linear3DFiniteElement::Linear3DFiniteElement()
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 
1230 void Linear3DFiniteElement::CalcShape(const IntegrationPoint &ip,
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 
1239 void Linear3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
1240  DenseMatrix &dshape) const
1241 {
1242  if (dshape.Height() == 4)
1243  {
1244  double *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 
1259 void Linear3DFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
1260 const
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.
1270 LinearWedgeFiniteElement::LinearWedgeFiniteElement()
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 
1293 void LinearWedgeFiniteElement::CalcShape(const IntegrationPoint &ip,
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 
1304 void LinearWedgeFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
1332 void LinearWedgeFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
1333 const
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 
1343 LinearPyramidFiniteElement::LinearPyramidFiniteElement()
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 
1363 void LinearPyramidFiniteElement::CalcShape(const IntegrationPoint &ip,
1364  Vector &shape) const
1365 {
1366  double x = ip.x, y = ip.y, z = ip.z;
1367  double ox = 1.-x-z, oy = 1.-y-z, oz = 1.-z;
1368 
1369  double 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  double 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 
1393 void LinearPyramidFiniteElement::CalcDShape(const IntegrationPoint &ip,
1394  DenseMatrix &dshape) const
1395 {
1396  double x = ip.x, y = ip.y, z = ip.z;
1397  double ox = 1.-x-z, oy = 1.-y-z, oz = 1.-z;
1398 
1399  double 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  double 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 
1454 void LinearPyramidFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
1455 const
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 
1465 Quadratic3DFiniteElement::Quadratic3DFiniteElement()
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 
1500 void Quadratic3DFiniteElement::CalcShape(const IntegrationPoint &ip,
1501  Vector &shape) const
1502 {
1503  double 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 
1522 void Quadratic3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
1523  DenseMatrix &dshape) const
1524 {
1525  double 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 
1544 TriLinear3DFiniteElement::TriLinear3DFiniteElement()
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 
1580 void TriLinear3DFiniteElement::CalcShape(const IntegrationPoint &ip,
1581  Vector &shape) const
1582 {
1583  double x = ip.x, y = ip.y, z = ip.z;
1584  double 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 
1596 void TriLinear3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
1597  DenseMatrix &dshape) const
1598 {
1599  double x = ip.x, y = ip.y, z = ip.z;
1600  double 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 
1636 P0SegmentFiniteElement::P0SegmentFiniteElement(int Ord)
1637  : NodalFiniteElement(1, Geometry::SEGMENT, 1, Ord) // default Ord = 0
1638 {
1639  Nodes.IntPoint(0).x = 0.5;
1640 }
1641 
1642 void P0SegmentFiniteElement::CalcShape(const IntegrationPoint &ip,
1643  Vector &shape) const
1644 {
1645  shape(0) = 1.0;
1646 }
1647 
1648 void P0SegmentFiniteElement::CalcDShape(const IntegrationPoint &ip,
1649  DenseMatrix &dshape) const
1650 {
1651  dshape(0,0) = 0.0;
1652 }
1653 
1654 CrouzeixRaviartFiniteElement::CrouzeixRaviartFiniteElement()
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 
1665 void CrouzeixRaviartFiniteElement::CalcShape(const IntegrationPoint &ip,
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 
1673 void CrouzeixRaviartFiniteElement::CalcDShape(const IntegrationPoint &ip,
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 
1681 CrouzeixRaviartQuadFiniteElement::CrouzeixRaviartQuadFiniteElement()
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 
1696 void CrouzeixRaviartQuadFiniteElement::CalcShape(const IntegrationPoint &ip,
1697  Vector &shape) const
1698 {
1699  const double 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 
1707 void CrouzeixRaviartQuadFiniteElement::CalcDShape(const IntegrationPoint &ip,
1708  DenseMatrix &dshape) const
1709 {
1710  const double 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 
1719 RT0TriangleFiniteElement::RT0TriangleFiniteElement()
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 
1730 void RT0TriangleFiniteElement::CalcVShape(const IntegrationPoint &ip,
1731  DenseMatrix &shape) const
1732 {
1733  double 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 
1743 void RT0TriangleFiniteElement::CalcDivShape(const IntegrationPoint &ip,
1744  Vector &divshape) const
1745 {
1746  divshape(0) = 2.;
1747  divshape(1) = 2.;
1748  divshape(2) = 2.;
1749 }
1750 
1751 const double RT0TriangleFiniteElement::nk[3][2] =
1752 { {0, -1}, {1, 1}, {-1, 0} };
1753 
1754 void RT0TriangleFiniteElement::GetLocalInterpolation (
1755  ElementTransformation &Trans, DenseMatrix &I) const
1756 {
1757  int k, j;
1758 #ifdef MFEM_THREAD_SAFE
1759  DenseMatrix vshape(dof, dim);
1760 #endif
1761 
1762 #ifdef MFEM_DEBUG
1763  for (k = 0; k < 3; k++)
1764  {
1765  CalcVShape (Nodes.IntPoint(k), vshape);
1766  for (j = 0; j < 3; j++)
1767  {
1768  double 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 
1780  IntegrationPoint ip;
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  double 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 
1806 void RT0TriangleFiniteElement::Project (
1807  VectorCoefficient &vc, ElementTransformation &Trans,
1808  Vector &dofs) const
1809 {
1810  double 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 
1826 RT0QuadFiniteElement::RT0QuadFiniteElement()
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 
1840 void RT0QuadFiniteElement::CalcVShape(const IntegrationPoint &ip,
1841  DenseMatrix &shape) const
1842 {
1843  double 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 
1855 void RT0QuadFiniteElement::CalcDivShape(const IntegrationPoint &ip,
1856  Vector &divshape) const
1857 {
1858  divshape(0) = 1.;
1859  divshape(1) = 1.;
1860  divshape(2) = 1.;
1861  divshape(3) = 1.;
1862 }
1863 
1864 const double RT0QuadFiniteElement::nk[4][2] =
1865 { {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
1866 
1867 void RT0QuadFiniteElement::GetLocalInterpolation (
1868  ElementTransformation &Trans, DenseMatrix &I) const
1869 {
1870  int k, j;
1871 #ifdef MFEM_THREAD_SAFE
1872  DenseMatrix vshape(dof, dim);
1873 #endif
1874 
1875 #ifdef MFEM_DEBUG
1876  for (k = 0; k < 4; k++)
1877  {
1878  CalcVShape (Nodes.IntPoint(k), vshape);
1879  for (j = 0; j < 4; j++)
1880  {
1881  double 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  {
1885  mfem::err << "RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
1886  " k = " << k << ", j = " << j << ", d = " << d << endl;
1887  mfem_error();
1888  }
1889  }
1890  }
1891 #endif
1892 
1893  IntegrationPoint ip;
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  double 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 
1919 void RT0QuadFiniteElement::Project (
1920  VectorCoefficient &vc, ElementTransformation &Trans,
1921  Vector &dofs) const
1922 {
1923  double 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 
1939 RT1TriangleFiniteElement::RT1TriangleFiniteElement()
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 
1960 void RT1TriangleFiniteElement::CalcVShape(const IntegrationPoint &ip,
1961  DenseMatrix &shape) const
1962 {
1963  double 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 
1983 void RT1TriangleFiniteElement::CalcDivShape(const IntegrationPoint &ip,
1984  Vector &divshape) const
1985 {
1986  double 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 
1998 const double 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 
2006 void RT1TriangleFiniteElement::GetLocalInterpolation (
2007  ElementTransformation &Trans, DenseMatrix &I) const
2008 {
2009  int k, j;
2010 #ifdef MFEM_THREAD_SAFE
2011  DenseMatrix vshape(dof, dim);
2012 #endif
2013 
2014 #ifdef MFEM_DEBUG
2015  for (k = 0; k < 8; k++)
2016  {
2017  CalcVShape (Nodes.IntPoint(k), vshape);
2018  for (j = 0; j < 8; j++)
2019  {
2020  double 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  {
2024  mfem::err << "RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2025  " k = " << k << ", j = " << j << ", d = " << d << endl;
2026  mfem_error();
2027  }
2028  }
2029  }
2030 #endif
2031 
2032  IntegrationPoint ip;
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  double 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 
2058 void RT1TriangleFiniteElement::Project (
2059  VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
2060 {
2061  double 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 
2078 RT1QuadFiniteElement::RT1QuadFiniteElement()
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 
2114 void RT1QuadFiniteElement::CalcVShape(const IntegrationPoint &ip,
2115  DenseMatrix &shape) const
2116 {
2117  double 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 
2151 void RT1QuadFiniteElement::CalcDivShape(const IntegrationPoint &ip,
2152  Vector &divshape) const
2153 {
2154  double 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 
2170 const double RT1QuadFiniteElement::nk[12][2] =
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 
2186 void RT1QuadFiniteElement::GetLocalInterpolation (
2187  ElementTransformation &Trans, DenseMatrix &I) const
2188 {
2189  int k, j;
2190 #ifdef MFEM_THREAD_SAFE
2191  DenseMatrix vshape(dof, dim);
2192 #endif
2193 
2194 #ifdef MFEM_DEBUG
2195  for (k = 0; k < 12; k++)
2196  {
2197  CalcVShape (Nodes.IntPoint(k), vshape);
2198  for (j = 0; j < 12; j++)
2199  {
2200  double 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  {
2204  mfem::err << "RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
2205  " k = " << k << ", j = " << j << ", d = " << d << endl;
2206  mfem_error();
2207  }
2208  }
2209  }
2210 #endif
2211 
2212  IntegrationPoint ip;
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  double 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 
2238 void RT1QuadFiniteElement::Project (
2239  VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
2240 {
2241  double 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 
2257 const double 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 
2316 RT2TriangleFiniteElement::RT2TriangleFiniteElement()
2317  : VectorFiniteElement(2, Geometry::TRIANGLE, 15, 3, H_DIV)
2318 {
2319  const double 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 
2353 void RT2TriangleFiniteElement::CalcVShape(const IntegrationPoint &ip,
2354  DenseMatrix &shape) const
2355 {
2356  double x = ip.x, y = ip.y;
2357 
2358  double 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  double 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  double 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 
2378 void RT2TriangleFiniteElement::CalcDivShape(const IntegrationPoint &ip,
2379  Vector &divshape) const
2380 {
2381  double x = ip.x, y = ip.y;
2382 
2383  double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
2384  4.*x*x, 4.*x*y, 4.*y*y
2385  };
2386 
2387  for (int i = 0; i < 15; i++)
2388  {
2389  double div = 0.0;
2390  for (int j = 0; j < 15; j++)
2391  {
2392  div += M[i][j] * DivB[j];
2393  }
2394  divshape(i) = div;
2395  }
2396 }
2397 
2398 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
2399 
2400 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
2401 
2402 RT2QuadFiniteElement::RT2QuadFiniteElement()
2403  : VectorFiniteElement(2, Geometry::SQUARE, 24, 3, H_DIV,
2404  FunctionSpace::Qk)
2405 {
2406  // y = 0 (pt[0])
2407  Nodes.IntPoint(0).x = dpt[0]; Nodes.IntPoint(0).y = pt[0];
2408  Nodes.IntPoint(1).x = dpt[1]; Nodes.IntPoint(1).y = pt[0];
2409  Nodes.IntPoint(2).x = dpt[2]; Nodes.IntPoint(2).y = pt[0];
2410  // x = 1 (pt[3])
2411  Nodes.IntPoint(3).x = pt[3]; Nodes.IntPoint(3).y = dpt[0];
2412  Nodes.IntPoint(4).x = pt[3]; Nodes.IntPoint(4).y = dpt[1];
2413  Nodes.IntPoint(5).x = pt[3]; Nodes.IntPoint(5).y = dpt[2];
2414  // y = 1 (pt[3])
2415  Nodes.IntPoint(6).x = dpt[2]; Nodes.IntPoint(6).y = pt[3];
2416  Nodes.IntPoint(7).x = dpt[1]; Nodes.IntPoint(7).y = pt[3];
2417  Nodes.IntPoint(8).x = dpt[0]; Nodes.IntPoint(8).y = pt[3];
2418  // x = 0 (pt[0])
2419  Nodes.IntPoint(9).x = pt[0]; Nodes.IntPoint(9).y = dpt[2];
2420  Nodes.IntPoint(10).x = pt[0]; Nodes.IntPoint(10).y = dpt[1];
2421  Nodes.IntPoint(11).x = pt[0]; Nodes.IntPoint(11).y = dpt[0];
2422  // x = pt[1] (interior)
2423  Nodes.IntPoint(12).x = pt[1]; Nodes.IntPoint(12).y = dpt[0];
2424  Nodes.IntPoint(13).x = pt[1]; Nodes.IntPoint(13).y = dpt[1];
2425  Nodes.IntPoint(14).x = pt[1]; Nodes.IntPoint(14).y = dpt[2];
2426  // x = pt[2] (interior)
2427  Nodes.IntPoint(15).x = pt[2]; Nodes.IntPoint(15).y = dpt[0];
2428  Nodes.IntPoint(16).x = pt[2]; Nodes.IntPoint(16).y = dpt[1];
2429  Nodes.IntPoint(17).x = pt[2]; Nodes.IntPoint(17).y = dpt[2];
2430  // y = pt[1] (interior)
2431  Nodes.IntPoint(18).x = dpt[0]; Nodes.IntPoint(18).y = pt[1];
2432  Nodes.IntPoint(19).x = dpt[1]; Nodes.IntPoint(19).y = pt[1];
2433  Nodes.IntPoint(20).x = dpt[2]; Nodes.IntPoint(20).y = pt[1];
2434  // y = pt[2] (interior)
2435  Nodes.IntPoint(21).x = dpt[0]; Nodes.IntPoint(21).y = pt[2];
2436  Nodes.IntPoint(22).x = dpt[1]; Nodes.IntPoint(22).y = pt[2];
2437  Nodes.IntPoint(23).x = dpt[2]; Nodes.IntPoint(23).y = pt[2];
2438 }
2439 
2440 void RT2QuadFiniteElement::CalcVShape(const IntegrationPoint &ip,
2441  DenseMatrix &shape) const
2442 {
2443  double x = ip.x, y = ip.y;
2444 
2445  double ax0 = pt[0] - x;
2446  double ax1 = pt[1] - x;
2447  double ax2 = pt[2] - x;
2448  double ax3 = pt[3] - x;
2449 
2450  double by0 = dpt[0] - y;
2451  double by1 = dpt[1] - y;
2452  double by2 = dpt[2] - y;
2453 
2454  double ay0 = pt[0] - y;
2455  double ay1 = pt[1] - y;
2456  double ay2 = pt[2] - y;
2457  double ay3 = pt[3] - y;
2458 
2459  double bx0 = dpt[0] - x;
2460  double bx1 = dpt[1] - x;
2461  double bx2 = dpt[2] - x;
2462 
2463  double A01 = pt[0] - pt[1];
2464  double A02 = pt[0] - pt[2];
2465  double A12 = pt[1] - pt[2];
2466  double A03 = pt[0] - pt[3];
2467  double A13 = pt[1] - pt[3];
2468  double A23 = pt[2] - pt[3];
2469 
2470  double B01 = dpt[0] - dpt[1];
2471  double B02 = dpt[0] - dpt[2];
2472  double B12 = dpt[1] - dpt[2];
2473 
2474  double tx0 = (bx1*bx2)/(B01*B02);
2475  double tx1 = -(bx0*bx2)/(B01*B12);
2476  double tx2 = (bx0*bx1)/(B02*B12);
2477 
2478  double ty0 = (by1*by2)/(B01*B02);
2479  double ty1 = -(by0*by2)/(B01*B12);
2480  double ty2 = (by0*by1)/(B02*B12);
2481 
2482  // y = 0 (p[0])
2483  shape(0, 0) = 0;
2484  shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
2485  shape(1, 0) = 0;
2486  shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
2487  shape(2, 0) = 0;
2488  shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
2489  // x = 1 (p[3])
2490  shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
2491  shape(3, 1) = 0;
2492  shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
2493  shape(4, 1) = 0;
2494  shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
2495  shape(5, 1) = 0;
2496  // y = 1 (p[3])
2497  shape(6, 0) = 0;
2498  shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
2499  shape(7, 0) = 0;
2500  shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
2501  shape(8, 0) = 0;
2502  shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
2503  // x = 0 (p[0])
2504  shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
2505  shape(9, 1) = 0;
2506  shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
2507  shape(10, 1) = 0;
2508  shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
2509  shape(11, 1) = 0;
2510  // x = p[1] (interior)
2511  shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
2512  shape(12, 1) = 0;
2513  shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
2514  shape(13, 1) = 0;
2515  shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
2516  shape(14, 1) = 0;
2517  // x = p[2] (interior)
2518  shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
2519  shape(15, 1) = 0;
2520  shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
2521  shape(16, 1) = 0;
2522  shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
2523  shape(17, 1) = 0;
2524  // y = p[1] (interior)
2525  shape(18, 0) = 0;
2526  shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
2527  shape(19, 0) = 0;
2528  shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
2529  shape(20, 0) = 0;
2530  shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
2531  // y = p[2] (interior)
2532  shape(21, 0) = 0;
2533  shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
2534  shape(22, 0) = 0;
2535  shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
2536  shape(23, 0) = 0;
2537  shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
2538 }
2539 
2540 void RT2QuadFiniteElement::CalcDivShape(const IntegrationPoint &ip,
2541  Vector &divshape) const
2542 {
2543  double x = ip.x, y = ip.y;
2544 
2545  double a01 = pt[0]*pt[1];
2546  double a02 = pt[0]*pt[2];
2547  double a12 = pt[1]*pt[2];
2548  double a03 = pt[0]*pt[3];
2549  double a13 = pt[1]*pt[3];
2550  double a23 = pt[2]*pt[3];
2551 
2552  double bx0 = dpt[0] - x;
2553  double bx1 = dpt[1] - x;
2554  double bx2 = dpt[2] - x;
2555 
2556  double by0 = dpt[0] - y;
2557  double by1 = dpt[1] - y;
2558  double by2 = dpt[2] - y;
2559 
2560  double A01 = pt[0] - pt[1];
2561  double A02 = pt[0] - pt[2];
2562  double A12 = pt[1] - pt[2];
2563  double A03 = pt[0] - pt[3];
2564  double A13 = pt[1] - pt[3];
2565  double A23 = pt[2] - pt[3];
2566 
2567  double A012 = pt[0] + pt[1] + pt[2];
2568  double A013 = pt[0] + pt[1] + pt[3];
2569  double A023 = pt[0] + pt[2] + pt[3];
2570  double A123 = pt[1] + pt[2] + pt[3];
2571 
2572  double B01 = dpt[0] - dpt[1];
2573  double B02 = dpt[0] - dpt[2];
2574  double B12 = dpt[1] - dpt[2];
2575 
2576  double tx0 = (bx1*bx2)/(B01*B02);
2577  double tx1 = -(bx0*bx2)/(B01*B12);
2578  double tx2 = (bx0*bx1)/(B02*B12);
2579 
2580  double ty0 = (by1*by2)/(B01*B02);
2581  double ty1 = -(by0*by2)/(B01*B12);
2582  double ty2 = (by0*by1)/(B02*B12);
2583 
2584  // y = 0 (p[0])
2585  divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
2586  divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
2587  divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
2588  // x = 1 (p[3])
2589  divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
2590  divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
2591  divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
2592  // y = 1 (p[3])
2593  divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
2594  divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
2595  divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
2596  // x = 0 (p[0])
2597  divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
2598  divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
2599  divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
2600  // x = p[1] (interior)
2601  divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
2602  divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
2603  divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
2604  // x = p[2] (interior)
2605  divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
2606  divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
2607  divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
2608  // y = p[1] (interior)
2609  divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
2610  divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
2611  divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
2612  // y = p[2] (interior)
2613  divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
2614  divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
2615  divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
2616 }
2617 
2618 const double RT2QuadFiniteElement::nk[24][2] =
2619 {
2620  // y = 0
2621  {0,-1}, {0,-1}, {0,-1},
2622  // x = 1
2623  {1, 0}, {1, 0}, {1, 0},
2624  // y = 1
2625  {0, 1}, {0, 1}, {0, 1},
2626  // x = 0
2627  {-1,0}, {-1,0}, {-1,0},
2628  // x = p[1] (interior)
2629  {1, 0}, {1, 0}, {1, 0},
2630  // x = p[2] (interior)
2631  {1, 0}, {1, 0}, {1, 0},
2632  // y = p[1] (interior)
2633  {0, 1}, {0, 1}, {0, 1},
2634  // y = p[1] (interior)
2635  {0, 1}, {0, 1}, {0, 1}
2636 };
2637 
2638 void RT2QuadFiniteElement::GetLocalInterpolation (
2639  ElementTransformation &Trans, DenseMatrix &I) const
2640 {
2641  int k, j;
2642 #ifdef MFEM_THREAD_SAFE
2643  DenseMatrix vshape(dof, dim);
2644 #endif
2645 
2646 #ifdef MFEM_DEBUG
2647  for (k = 0; k < 24; k++)
2648  {
2649  CalcVShape (Nodes.IntPoint(k), vshape);
2650  for (j = 0; j < 24; j++)
2651  {
2652  double d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
2653  if (j == k) { d -= 1.0; }
2654  if (fabs(d) > 1.0e-12)
2655  {
2656  mfem::err << "RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
2657  " k = " << k << ", j = " << j << ", d = " << d << endl;
2658  mfem_error();
2659  }
2660  }
2661  }
2662 #endif
2663 
2664  IntegrationPoint ip;
2665  ip.x = ip.y = 0.0;
2666  Trans.SetIntPoint (&ip);
2667  // Trans must be linear (more to have embedding?)
2668  // set Jinv = |J| J^{-t} = adj(J)^t
2669  const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2670 
2671  double vk[2];
2672  Vector xk (vk, 2);
2673 
2674  for (k = 0; k < 24; k++)
2675  {
2676  Trans.Transform (Nodes.IntPoint (k), xk);
2677  ip.x = vk[0]; ip.y = vk[1];
2678  CalcVShape (ip, vshape);
2679  // vk = |J| J^{-t} nk
2680  vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
2681  vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
2682  for (j = 0; j < 24; j++)
2683  if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
2684  {
2685  I(k,j) = 0.0;
2686  }
2687  }
2688 }
2689 
2690 void RT2QuadFiniteElement::Project (
2691  VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
2692 {
2693  double vk[2];
2694  Vector xk (vk, 2);
2695 
2696  for (int k = 0; k < 24; k++)
2697  {
2698  Trans.SetIntPoint (&Nodes.IntPoint (k));
2699  // set Jinv = |J| J^{-t} = adj(J)^t
2700  const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
2701 
2702  vc.Eval (xk, Trans, Nodes.IntPoint (k));
2703  // xk^t |J| J^{-t} nk
2704  dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
2705  vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
2706  }
2707 }
2708 
2709 P1SegmentFiniteElement::P1SegmentFiniteElement()
2710  : NodalFiniteElement(1, Geometry::SEGMENT, 2, 1)
2711 {
2712  Nodes.IntPoint(0).x = 0.33333333333333333333;
2713  Nodes.IntPoint(1).x = 0.66666666666666666667;
2714 }
2715 
2716 void P1SegmentFiniteElement::CalcShape(const IntegrationPoint &ip,
2717  Vector &shape) const
2718 {
2719  double x = ip.x;
2720 
2721  shape(0) = 2. - 3. * x;
2722  shape(1) = 3. * x - 1.;
2723 }
2724 
2725 void P1SegmentFiniteElement::CalcDShape(const IntegrationPoint &ip,
2726  DenseMatrix &dshape) const
2727 {
2728  dshape(0,0) = -3.;
2729  dshape(1,0) = 3.;
2730 }
2731 
2732 
2733 P2SegmentFiniteElement::P2SegmentFiniteElement()
2734  : NodalFiniteElement(1, Geometry::SEGMENT, 3, 2)
2735 {
2736  const double p = 0.11270166537925831148;
2737 
2738  Nodes.IntPoint(0).x = p;
2739  Nodes.IntPoint(1).x = 0.5;
2740  Nodes.IntPoint(2).x = 1.-p;
2741 }
2742 
2743 void P2SegmentFiniteElement::CalcShape(const IntegrationPoint &ip,
2744  Vector &shape) const
2745 {
2746  const double p = 0.11270166537925831148;
2747  const double w = 1./((1-2*p)*(1-2*p));
2748  double x = ip.x;
2749 
2750  shape(0) = (2*x-1)*(x-1+p)*w;
2751  shape(1) = 4*(x-1+p)*(p-x)*w;
2752  shape(2) = (2*x-1)*(x-p)*w;
2753 }
2754 
2755 void P2SegmentFiniteElement::CalcDShape(const IntegrationPoint &ip,
2756  DenseMatrix &dshape) const
2757 {
2758  const double p = 0.11270166537925831148;
2759  const double w = 1./((1-2*p)*(1-2*p));
2760  double x = ip.x;
2761 
2762  dshape(0,0) = (-3+4*x+2*p)*w;
2763  dshape(1,0) = (4-8*x)*w;
2764  dshape(2,0) = (-1+4*x-2*p)*w;
2765 }
2766 
2767 
2768 Lagrange1DFiniteElement::Lagrange1DFiniteElement(int degree)
2769  : NodalFiniteElement(1, Geometry::SEGMENT, degree+1, degree)
2770 {
2771  int i, m = degree;
2772 
2773  Nodes.IntPoint(0).x = 0.0;
2774  Nodes.IntPoint(1).x = 1.0;
2775  for (i = 1; i < m; i++)
2776  {
2777  Nodes.IntPoint(i+1).x = double(i) / m;
2778  }
2779 
2780  rwk.SetSize(degree+1);
2781 #ifndef MFEM_THREAD_SAFE
2782  rxxk.SetSize(degree+1);
2783 #endif
2784 
2785  rwk(0) = 1.0;
2786  for (i = 1; i <= m; i++)
2787  {
2788  rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
2789  }
2790  for (i = 0; i < m/2+1; i++)
2791  {
2792  rwk(m-i) = ( rwk(i) *= rwk(m-i) );
2793  }
2794  for (i = m-1; i >= 0; i -= 2)
2795  {
2796  rwk(i) = -rwk(i);
2797  }
2798 }
2799 
2800 void Lagrange1DFiniteElement::CalcShape(const IntegrationPoint &ip,
2801  Vector &shape) const
2802 {
2803  double w, wk, x = ip.x;
2804  int i, k, m = GetOrder();
2805 
2806 #ifdef MFEM_THREAD_SAFE
2807  Vector rxxk(m+1);
2808 #endif
2809 
2810  k = (int) floor ( m * x + 0.5 );
2811  k = k > m ? m : k < 0 ? 0 : k; // clamp k to [0,m]
2812 
2813  wk = 1.0;
2814  for (i = 0; i <= m; i++)
2815  if (i != k)
2816  {
2817  wk *= ( rxxk(i) = x - (double)(i) / m );
2818  }
2819  w = wk * ( rxxk(k) = x - (double)(k) / m );
2820 
2821  if (k != 0)
2822  {
2823  shape(0) = w * rwk(0) / rxxk(0);
2824  }
2825  else
2826  {
2827  shape(0) = wk * rwk(0);
2828  }
2829  if (k != m)
2830  {
2831  shape(1) = w * rwk(m) / rxxk(m);
2832  }
2833  else
2834  {
2835  shape(1) = wk * rwk(k);
2836  }
2837  for (i = 1; i < m; i++)
2838  if (i != k)
2839  {
2840  shape(i+1) = w * rwk(i) / rxxk(i);
2841  }
2842  else
2843  {
2844  shape(k+1) = wk * rwk(k);
2845  }
2846 }
2847 
2848 void Lagrange1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
2849  DenseMatrix &dshape) const
2850 {
2851  double s, srx, w, wk, x = ip.x;
2852  int i, k, m = GetOrder();
2853 
2854 #ifdef MFEM_THREAD_SAFE
2855  Vector rxxk(m+1);
2856 #endif
2857 
2858  k = (int) floor ( m * x + 0.5 );
2859  k = k > m ? m : k < 0 ? 0 : k; // clamp k to [0,m]
2860 
2861  wk = 1.0;
2862  for (i = 0; i <= m; i++)
2863  if (i != k)
2864  {
2865  wk *= ( rxxk(i) = x - (double)(i) / m );
2866  }
2867  w = wk * ( rxxk(k) = x - (double)(k) / m );
2868 
2869  for (i = 0; i <= m; i++)
2870  {
2871  rxxk(i) = 1.0 / rxxk(i);
2872  }
2873  srx = 0.0;
2874  for (i = 0; i <= m; i++)
2875  if (i != k)
2876  {
2877  srx += rxxk(i);
2878  }
2879  s = w * srx + wk;
2880 
2881  if (k != 0)
2882  {
2883  dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
2884  }
2885  else
2886  {
2887  dshape(0,0) = wk * srx * rwk(0);
2888  }
2889  if (k != m)
2890  {
2891  dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
2892  }
2893  else
2894  {
2895  dshape(1,0) = wk * srx * rwk(k);
2896  }
2897  for (i = 1; i < m; i++)
2898  if (i != k)
2899  {
2900  dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
2901  }
2902  else
2903  {
2904  dshape(k+1,0) = wk * srx * rwk(k);
2905  }
2906 }
2907 
2908 
2909 P1TetNonConfFiniteElement::P1TetNonConfFiniteElement()
2910  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 4, 1)
2911 {
2912  Nodes.IntPoint(0).x = 0.33333333333333333333;
2913  Nodes.IntPoint(0).y = 0.33333333333333333333;
2914  Nodes.IntPoint(0).z = 0.33333333333333333333;
2915 
2916  Nodes.IntPoint(1).x = 0.0;
2917  Nodes.IntPoint(1).y = 0.33333333333333333333;
2918  Nodes.IntPoint(1).z = 0.33333333333333333333;
2919 
2920  Nodes.IntPoint(2).x = 0.33333333333333333333;
2921  Nodes.IntPoint(2).y = 0.0;
2922  Nodes.IntPoint(2).z = 0.33333333333333333333;
2923 
2924  Nodes.IntPoint(3).x = 0.33333333333333333333;
2925  Nodes.IntPoint(3).y = 0.33333333333333333333;
2926  Nodes.IntPoint(3).z = 0.0;
2927 
2928 }
2929 
2930 void P1TetNonConfFiniteElement::CalcShape(const IntegrationPoint &ip,
2931  Vector &shape) const
2932 {
2933  double L0, L1, L2, L3;
2934 
2935  L1 = ip.x; L2 = ip.y; L3 = ip.z; L0 = 1.0 - L1 - L2 - L3;
2936  shape(0) = 1.0 - 3.0 * L0;
2937  shape(1) = 1.0 - 3.0 * L1;
2938  shape(2) = 1.0 - 3.0 * L2;
2939  shape(3) = 1.0 - 3.0 * L3;
2940 }
2941 
2942 void P1TetNonConfFiniteElement::CalcDShape(const IntegrationPoint &ip,
2943  DenseMatrix &dshape) const
2944 {
2945  dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
2946  dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
2947  dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
2948  dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
2949 }
2950 
2951 
2952 P0TetFiniteElement::P0TetFiniteElement()
2953  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 1, 0)
2954 {
2955  Nodes.IntPoint(0).x = 0.25;
2956  Nodes.IntPoint(0).y = 0.25;
2957  Nodes.IntPoint(0).z = 0.25;
2958 }
2959 
2960 void P0TetFiniteElement::CalcShape(const IntegrationPoint &ip,
2961  Vector &shape) const
2962 {
2963  shape(0) = 1.0;
2964 }
2965 
2966 void P0TetFiniteElement::CalcDShape(const IntegrationPoint &ip,
2967  DenseMatrix &dshape) const
2968 {
2969  dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
2970 }
2971 
2972 
2973 P0HexFiniteElement::P0HexFiniteElement()
2974  : NodalFiniteElement(3, Geometry::CUBE, 1, 0, FunctionSpace::Qk)
2975 {
2976  Nodes.IntPoint(0).x = 0.5;
2977  Nodes.IntPoint(0).y = 0.5;
2978  Nodes.IntPoint(0).z = 0.5;
2979 }
2980 
2981 void P0HexFiniteElement::CalcShape(const IntegrationPoint &ip,
2982  Vector &shape) const
2983 {
2984  shape(0) = 1.0;
2985 }
2986 
2987 void P0HexFiniteElement::CalcDShape(const IntegrationPoint &ip,
2988  DenseMatrix &dshape) const
2989 {
2990  dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
2991 }
2992 
2993 
2994 P0WdgFiniteElement::P0WdgFiniteElement()
2995  : NodalFiniteElement(3, Geometry::PRISM, 1, 0, FunctionSpace::Qk)
2996 {
2997  Nodes.IntPoint(0).x = 1.0 / 3.0;
2998  Nodes.IntPoint(0).y = 1.0 / 3.0;
2999  Nodes.IntPoint(0).z = 0.5;
3000 }
3001 
3002 void P0WdgFiniteElement::CalcShape(const IntegrationPoint &ip,
3003  Vector &shape) const
3004 {
3005  shape(0) = 1.0;
3006 }
3007 
3008 void P0WdgFiniteElement::CalcDShape(const IntegrationPoint &ip,
3009  DenseMatrix &dshape) const
3010 {
3011  dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3012 }
3013 
3014 
3015 P0PyrFiniteElement::P0PyrFiniteElement()
3016  : NodalFiniteElement(3, Geometry::PYRAMID, 1, 0, FunctionSpace::Qk)
3017 {
3018  Nodes.IntPoint(0).x = 0.375;
3019  Nodes.IntPoint(0).y = 0.375;
3020  Nodes.IntPoint(0).z = 0.25;
3021 }
3022 
3023 void P0PyrFiniteElement::CalcShape(const IntegrationPoint &ip,
3024  Vector &shape) const
3025 {
3026  shape(0) = 1.0;
3027 }
3028 
3029 void P0PyrFiniteElement::CalcDShape(const IntegrationPoint &ip,
3030  DenseMatrix &dshape) const
3031 {
3032  dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
3033 }
3034 
3035 
3036 LagrangeHexFiniteElement::LagrangeHexFiniteElement (int degree)
3037  : NodalFiniteElement(3, Geometry::CUBE, (degree+1)*(degree+1)*(degree+1),
3038  degree, FunctionSpace::Qk)
3039 {
3040  if (degree == 2)
3041  {
3042  I = new int[dof];
3043  J = new int[dof];
3044  K = new int[dof];
3045  // nodes
3046  I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3047  I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3048  I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3049  I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3050  I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3051  I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3052  I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3053  I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3054  // edges
3055  I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3056  I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
3057  I[10] = 2; J[10] = 1; K[10] = 0;
3058  I[11] = 0; J[11] = 2; K[11] = 0;
3059  I[12] = 2; J[12] = 0; K[12] = 1;
3060  I[13] = 1; J[13] = 2; K[13] = 1;
3061  I[14] = 2; J[14] = 1; K[14] = 1;
3062  I[15] = 0; J[15] = 2; K[15] = 1;
3063  I[16] = 0; J[16] = 0; K[16] = 2;
3064  I[17] = 1; J[17] = 0; K[17] = 2;
3065  I[18] = 1; J[18] = 1; K[18] = 2;
3066  I[19] = 0; J[19] = 1; K[19] = 2;
3067  // faces
3068  I[20] = 2; J[20] = 2; K[20] = 0;
3069  I[21] = 2; J[21] = 0; K[21] = 2;
3070  I[22] = 1; J[22] = 2; K[22] = 2;
3071  I[23] = 2; J[23] = 1; K[23] = 2;
3072  I[24] = 0; J[24] = 2; K[24] = 2;
3073  I[25] = 2; J[25] = 2; K[25] = 1;
3074  // element
3075  I[26] = 2; J[26] = 2; K[26] = 2;
3076  }
3077  else if (degree == 3)
3078  {
3079  I = new int[dof];
3080  J = new int[dof];
3081  K = new int[dof];
3082  // nodes
3083  I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
3084  I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
3085  I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
3086  I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
3087  I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
3088  I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
3089  I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
3090  I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
3091  // edges
3092  I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
3093  I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
3094  I[10] = 1; J[10] = 2; K[10] = 0;
3095  I[11] = 1; J[11] = 3; K[11] = 0;
3096  I[12] = 2; J[12] = 1; K[12] = 0;
3097  I[13] = 3; J[13] = 1; K[13] = 0;
3098  I[14] = 0; J[14] = 2; K[14] = 0;
3099  I[15] = 0; J[15] = 3; K[15] = 0;
3100  I[16] = 2; J[16] = 0; K[16] = 1;
3101  I[17] = 3; J[17] = 0; K[17] = 1;
3102  I[18] = 1; J[18] = 2; K[18] = 1;
3103  I[19] = 1; J[19] = 3; K[19] = 1;
3104  I[20] = 2; J[20] = 1; K[20] = 1;
3105  I[21] = 3; J[21] = 1; K[21] = 1;
3106  I[22] = 0; J[22] = 2; K[22] = 1;
3107  I[23] = 0; J[23] = 3; K[23] = 1;
3108  I[24] = 0; J[24] = 0; K[24] = 2;
3109  I[25] = 0; J[25] = 0; K[25] = 3;
3110  I[26] = 1; J[26] = 0; K[26] = 2;
3111  I[27] = 1; J[27] = 0; K[27] = 3;
3112  I[28] = 1; J[28] = 1; K[28] = 2;
3113  I[29] = 1; J[29] = 1; K[29] = 3;
3114  I[30] = 0; J[30] = 1; K[30] = 2;
3115  I[31] = 0; J[31] = 1; K[31] = 3;
3116  // faces
3117  I[32] = 2; J[32] = 3; K[32] = 0;
3118  I[33] = 3; J[33] = 3; K[33] = 0;
3119  I[34] = 2; J[34] = 2; K[34] = 0;
3120  I[35] = 3; J[35] = 2; K[35] = 0;
3121  I[36] = 2; J[36] = 0; K[36] = 2;
3122  I[37] = 3; J[37] = 0; K[37] = 2;
3123  I[38] = 2; J[38] = 0; K[38] = 3;
3124  I[39] = 3; J[39] = 0; K[39] = 3;
3125  I[40] = 1; J[40] = 2; K[40] = 2;
3126  I[41] = 1; J[41] = 3; K[41] = 2;
3127  I[42] = 1; J[42] = 2; K[42] = 3;
3128  I[43] = 1; J[43] = 3; K[43] = 3;
3129  I[44] = 3; J[44] = 1; K[44] = 2;
3130  I[45] = 2; J[45] = 1; K[45] = 2;
3131  I[46] = 3; J[46] = 1; K[46] = 3;
3132  I[47] = 2; J[47] = 1; K[47] = 3;
3133  I[48] = 0; J[48] = 3; K[48] = 2;
3134  I[49] = 0; J[49] = 2; K[49] = 2;
3135  I[50] = 0; J[50] = 3; K[50] = 3;
3136  I[51] = 0; J[51] = 2; K[51] = 3;
3137  I[52] = 2; J[52] = 2; K[52] = 1;
3138  I[53] = 3; J[53] = 2; K[53] = 1;
3139  I[54] = 2; J[54] = 3; K[54] = 1;
3140  I[55] = 3; J[55] = 3; K[55] = 1;
3141  // element
3142  I[56] = 2; J[56] = 2; K[56] = 2;
3143  I[57] = 3; J[57] = 2; K[57] = 2;
3144  I[58] = 3; J[58] = 3; K[58] = 2;
3145  I[59] = 2; J[59] = 3; K[59] = 2;
3146  I[60] = 2; J[60] = 2; K[60] = 3;
3147  I[61] = 3; J[61] = 2; K[61] = 3;
3148  I[62] = 3; J[62] = 3; K[62] = 3;
3149  I[63] = 2; J[63] = 3; K[63] = 3;
3150  }
3151  else
3152  {
3153  mfem_error ("LagrangeHexFiniteElement::LagrangeHexFiniteElement");
3154  }
3155 
3156  fe1d = new Lagrange1DFiniteElement(degree);
3157  dof1d = fe1d -> GetDof();
3158 
3159 #ifndef MFEM_THREAD_SAFE
3160  shape1dx.SetSize(dof1d);
3161  shape1dy.SetSize(dof1d);
3162  shape1dz.SetSize(dof1d);
3163 
3164  dshape1dx.SetSize(dof1d,1);
3165  dshape1dy.SetSize(dof1d,1);
3166  dshape1dz.SetSize(dof1d,1);
3167 #endif
3168 
3169  for (int n = 0; n < dof; n++)
3170  {
3171  Nodes.IntPoint(n).x = fe1d -> GetNodes().IntPoint(I[n]).x;
3172  Nodes.IntPoint(n).y = fe1d -> GetNodes().IntPoint(J[n]).x;
3173  Nodes.IntPoint(n).z = fe1d -> GetNodes().IntPoint(K[n]).x;
3174  }
3175 }
3176 
3177 void LagrangeHexFiniteElement::CalcShape(const IntegrationPoint &ip,
3178  Vector &shape) const
3179 {
3180  IntegrationPoint ipy, ipz;
3181  ipy.x = ip.y;
3182  ipz.x = ip.z;
3183 
3184 #ifdef MFEM_THREAD_SAFE
3185  Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3186 #endif
3187 
3188  fe1d -> CalcShape(ip, shape1dx);
3189  fe1d -> CalcShape(ipy, shape1dy);
3190  fe1d -> CalcShape(ipz, shape1dz);
3191 
3192  for (int n = 0; n < dof; n++)
3193  {
3194  shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
3195  }
3196 }
3197 
3198 void LagrangeHexFiniteElement::CalcDShape(const IntegrationPoint &ip,
3199  DenseMatrix &dshape) const
3200 {
3201  IntegrationPoint ipy, ipz;
3202  ipy.x = ip.y;
3203  ipz.x = ip.z;
3204 
3205 #ifdef MFEM_THREAD_SAFE
3206  Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
3207  DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
3208 #endif
3209 
3210  fe1d -> CalcShape(ip, shape1dx);
3211  fe1d -> CalcShape(ipy, shape1dy);
3212  fe1d -> CalcShape(ipz, shape1dz);
3213 
3214  fe1d -> CalcDShape(ip, dshape1dx);
3215  fe1d -> CalcDShape(ipy, dshape1dy);
3216  fe1d -> CalcDShape(ipz, dshape1dz);
3217 
3218  for (int n = 0; n < dof; n++)
3219  {
3220  dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
3221  dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
3222  dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
3223  }
3224 }
3225 
3226 LagrangeHexFiniteElement::~LagrangeHexFiniteElement ()
3227 {
3228  delete fe1d;
3229 
3230  delete [] I;
3231  delete [] J;
3232  delete [] K;
3233 }
3234 
3235 
3236 RefinedLinear1DFiniteElement::RefinedLinear1DFiniteElement()
3237  : NodalFiniteElement(1, Geometry::SEGMENT, 3, 4)
3238 {
3239  Nodes.IntPoint(0).x = 0.0;
3240  Nodes.IntPoint(1).x = 1.0;
3241  Nodes.IntPoint(2).x = 0.5;
3242 }
3243 
3244 void RefinedLinear1DFiniteElement::CalcShape(const IntegrationPoint &ip,
3245  Vector &shape) const
3246 {
3247  double x = ip.x;
3248 
3249  if (x <= 0.5)
3250  {
3251  shape(0) = 1.0 - 2.0 * x;
3252  shape(1) = 0.0;
3253  shape(2) = 2.0 * x;
3254  }
3255  else
3256  {
3257  shape(0) = 0.0;
3258  shape(1) = 2.0 * x - 1.0;
3259  shape(2) = 2.0 - 2.0 * x;
3260  }
3261 }
3262 
3263 void RefinedLinear1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
3264  DenseMatrix &dshape) const
3265 {
3266  double x = ip.x;
3267 
3268  if (x <= 0.5)
3269  {
3270  dshape(0,0) = - 2.0;
3271  dshape(1,0) = 0.0;
3272  dshape(2,0) = 2.0;
3273  }
3274  else
3275  {
3276  dshape(0,0) = 0.0;
3277  dshape(1,0) = 2.0;
3278  dshape(2,0) = - 2.0;
3279  }
3280 }
3281 
3282 RefinedLinear2DFiniteElement::RefinedLinear2DFiniteElement()
3283  : NodalFiniteElement(2, Geometry::TRIANGLE, 6, 5)
3284 {
3285  Nodes.IntPoint(0).x = 0.0;
3286  Nodes.IntPoint(0).y = 0.0;
3287  Nodes.IntPoint(1).x = 1.0;
3288  Nodes.IntPoint(1).y = 0.0;
3289  Nodes.IntPoint(2).x = 0.0;
3290  Nodes.IntPoint(2).y = 1.0;
3291  Nodes.IntPoint(3).x = 0.5;
3292  Nodes.IntPoint(3).y = 0.0;
3293  Nodes.IntPoint(4).x = 0.5;
3294  Nodes.IntPoint(4).y = 0.5;
3295  Nodes.IntPoint(5).x = 0.0;
3296  Nodes.IntPoint(5).y = 0.5;
3297 }
3298 
3299 void RefinedLinear2DFiniteElement::CalcShape(const IntegrationPoint &ip,
3300  Vector &shape) const
3301 {
3302  int i;
3303 
3304  double L0, L1, L2;
3305  L0 = 2.0 * ( 1. - ip.x - ip.y );
3306  L1 = 2.0 * ( ip.x );
3307  L2 = 2.0 * ( ip.y );
3308 
3309  // The reference triangle is split in 4 triangles as follows:
3310  //
3311  // T0 - 0,3,5
3312  // T1 - 1,3,4
3313  // T2 - 2,4,5
3314  // T3 - 3,4,5
3315 
3316  for (i = 0; i < 6; i++)
3317  {
3318  shape(i) = 0.0;
3319  }
3320 
3321  if (L0 >= 1.0) // T0
3322  {
3323  shape(0) = L0 - 1.0;
3324  shape(3) = L1;
3325  shape(5) = L2;
3326  }
3327  else if (L1 >= 1.0) // T1
3328  {
3329  shape(3) = L0;
3330  shape(1) = L1 - 1.0;
3331  shape(4) = L2;
3332  }
3333  else if (L2 >= 1.0) // T2
3334  {
3335  shape(5) = L0;
3336  shape(4) = L1;
3337  shape(2) = L2 - 1.0;
3338  }
3339  else // T3
3340  {
3341  shape(3) = 1.0 - L2;
3342  shape(4) = 1.0 - L0;
3343  shape(5) = 1.0 - L1;
3344  }
3345 }
3346 
3347 void RefinedLinear2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
3348  DenseMatrix &dshape) const
3349 {
3350  int i,j;
3351 
3352  double L0, L1, L2;
3353  L0 = 2.0 * ( 1. - ip.x - ip.y );
3354  L1 = 2.0 * ( ip.x );
3355  L2 = 2.0 * ( ip.y );
3356 
3357  double DL0[2], DL1[2], DL2[2];
3358  DL0[0] = -2.0; DL0[1] = -2.0;
3359  DL1[0] = 2.0; DL1[1] = 0.0;
3360  DL2[0] = 0.0; DL2[1] = 2.0;
3361 
3362  for (i = 0; i < 6; i++)
3363  for (j = 0; j < 2; j++)
3364  {
3365  dshape(i,j) = 0.0;
3366  }
3367 
3368  if (L0 >= 1.0) // T0
3369  {
3370  for (j = 0; j < 2; j++)
3371  {
3372  dshape(0,j) = DL0[j];
3373  dshape(3,j) = DL1[j];
3374  dshape(5,j) = DL2[j];
3375  }
3376  }
3377  else if (L1 >= 1.0) // T1
3378  {
3379  for (j = 0; j < 2; j++)
3380  {
3381  dshape(3,j) = DL0[j];
3382  dshape(1,j) = DL1[j];
3383  dshape(4,j) = DL2[j];
3384  }
3385  }
3386  else if (L2 >= 1.0) // T2
3387  {
3388  for (j = 0; j < 2; j++)
3389  {
3390  dshape(5,j) = DL0[j];
3391  dshape(4,j) = DL1[j];
3392  dshape(2,j) = DL2[j];
3393  }
3394  }
3395  else // T3
3396  {
3397  for (j = 0; j < 2; j++)
3398  {
3399  dshape(3,j) = - DL2[j];
3400  dshape(4,j) = - DL0[j];
3401  dshape(5,j) = - DL1[j];
3402  }
3403  }
3404 }
3405 
3406 RefinedLinear3DFiniteElement::RefinedLinear3DFiniteElement()
3407  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 10, 4)
3408 {
3409  Nodes.IntPoint(0).x = 0.0;
3410  Nodes.IntPoint(0).y = 0.0;
3411  Nodes.IntPoint(0).z = 0.0;
3412  Nodes.IntPoint(1).x = 1.0;
3413  Nodes.IntPoint(1).y = 0.0;
3414  Nodes.IntPoint(1).z = 0.0;
3415  Nodes.IntPoint(2).x = 0.0;
3416  Nodes.IntPoint(2).y = 1.0;
3417  Nodes.IntPoint(2).z = 0.0;
3418  Nodes.IntPoint(3).x = 0.0;
3419  Nodes.IntPoint(3).y = 0.0;
3420  Nodes.IntPoint(3).z = 1.0;
3421  Nodes.IntPoint(4).x = 0.5;
3422  Nodes.IntPoint(4).y = 0.0;
3423  Nodes.IntPoint(4).z = 0.0;
3424  Nodes.IntPoint(5).x = 0.0;
3425  Nodes.IntPoint(5).y = 0.5;
3426  Nodes.IntPoint(5).z = 0.0;
3427  Nodes.IntPoint(6).x = 0.0;
3428  Nodes.IntPoint(6).y = 0.0;
3429  Nodes.IntPoint(6).z = 0.5;
3430  Nodes.IntPoint(7).x = 0.5;
3431  Nodes.IntPoint(7).y = 0.5;
3432  Nodes.IntPoint(7).z = 0.0;
3433  Nodes.IntPoint(8).x = 0.5;
3434  Nodes.IntPoint(8).y = 0.0;
3435  Nodes.IntPoint(8).z = 0.5;
3436  Nodes.IntPoint(9).x = 0.0;
3437  Nodes.IntPoint(9).y = 0.5;
3438  Nodes.IntPoint(9).z = 0.5;
3439 }
3440 
3441 void RefinedLinear3DFiniteElement::CalcShape(const IntegrationPoint &ip,
3442  Vector &shape) const
3443 {
3444  int i;
3445 
3446  double L0, L1, L2, L3, L4, L5;
3447  L0 = 2.0 * ( 1. - ip.x - ip.y - ip.z );
3448  L1 = 2.0 * ( ip.x );
3449  L2 = 2.0 * ( ip.y );
3450  L3 = 2.0 * ( ip.z );
3451  L4 = 2.0 * ( ip.x + ip.y );
3452  L5 = 2.0 * ( ip.y + ip.z );
3453 
3454  // The reference tetrahedron is split in 8 tetrahedra as follows:
3455  //
3456  // T0 - 0,4,5,6
3457  // T1 - 1,4,7,8
3458  // T2 - 2,5,7,9
3459  // T3 - 3,6,8,9
3460  // T4 - 4,5,6,8
3461  // T5 - 4,5,7,8
3462  // T6 - 5,6,8,9
3463  // T7 - 5,7,8,9
3464 
3465  for (i = 0; i < 10; i++)
3466  {
3467  shape(i) = 0.0;
3468  }
3469 
3470  if (L0 >= 1.0) // T0
3471  {
3472  shape(0) = L0 - 1.0;
3473  shape(4) = L1;
3474  shape(5) = L2;
3475  shape(6) = L3;
3476  }
3477  else if (L1 >= 1.0) // T1
3478  {
3479  shape(4) = L0;
3480  shape(1) = L1 - 1.0;
3481  shape(7) = L2;
3482  shape(8) = L3;
3483  }
3484  else if (L2 >= 1.0) // T2
3485  {
3486  shape(5) = L0;
3487  shape(7) = L1;
3488  shape(2) = L2 - 1.0;
3489  shape(9) = L3;
3490  }
3491  else if (L3 >= 1.0) // T3
3492  {
3493  shape(6) = L0;
3494  shape(8) = L1;
3495  shape(9) = L2;
3496  shape(3) = L3 - 1.0;
3497  }
3498  else if ((L4 <= 1.0) && (L5 <= 1.0)) // T4
3499  {
3500  shape(4) = 1.0 - L5;
3501  shape(5) = L2;
3502  shape(6) = 1.0 - L4;
3503  shape(8) = 1.0 - L0;
3504  }
3505  else if ((L4 >= 1.0) && (L5 <= 1.0)) // T5
3506  {
3507  shape(4) = 1.0 - L5;
3508  shape(5) = 1.0 - L1;
3509  shape(7) = L4 - 1.0;
3510  shape(8) = L3;
3511  }
3512  else if ((L4 <= 1.0) && (L5 >= 1.0)) // T6
3513  {
3514  shape(5) = 1.0 - L3;
3515  shape(6) = 1.0 - L4;
3516  shape(8) = L1;
3517  shape(9) = L5 - 1.0;
3518  }
3519  else if ((L4 >= 1.0) && (L5 >= 1.0)) // T7
3520  {
3521  shape(5) = L0;
3522  shape(7) = L4 - 1.0;
3523  shape(8) = 1.0 - L2;
3524  shape(9) = L5 - 1.0;
3525  }
3526 }
3527 
3528 void RefinedLinear3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
3529  DenseMatrix &dshape) const
3530 {
3531  int i,j;
3532 
3533  double L0, L1, L2, L3, L4, L5;
3534  L0 = 2.0 * ( 1. - ip.x - ip.y - ip.z );
3535  L1 = 2.0 * ( ip.x );
3536  L2 = 2.0 * ( ip.y );
3537  L3 = 2.0 * ( ip.z );
3538  L4 = 2.0 * ( ip.x + ip.y );
3539  L5 = 2.0 * ( ip.y + ip.z );
3540 
3541  double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
3542  DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
3543  DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
3544  DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
3545  DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
3546  DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
3547  DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
3548 
3549  for (i = 0; i < 10; i++)
3550  for (j = 0; j < 3; j++)
3551  {
3552  dshape(i,j) = 0.0;
3553  }
3554 
3555  if (L0 >= 1.0) // T0
3556  {
3557  for (j = 0; j < 3; j++)
3558  {
3559  dshape(0,j) = DL0[j];
3560  dshape(4,j) = DL1[j];
3561  dshape(5,j) = DL2[j];
3562  dshape(6,j) = DL3[j];
3563  }
3564  }
3565  else if (L1 >= 1.0) // T1
3566  {
3567  for (j = 0; j < 3; j++)
3568  {
3569  dshape(4,j) = DL0[j];
3570  dshape(1,j) = DL1[j];
3571  dshape(7,j) = DL2[j];
3572  dshape(8,j) = DL3[j];
3573  }
3574  }
3575  else if (L2 >= 1.0) // T2
3576  {
3577  for (j = 0; j < 3; j++)
3578  {
3579  dshape(5,j) = DL0[j];
3580  dshape(7,j) = DL1[j];
3581  dshape(2,j) = DL2[j];
3582  dshape(9,j) = DL3[j];
3583  }
3584  }
3585  else if (L3 >= 1.0) // T3
3586  {
3587  for (j = 0; j < 3; j++)
3588  {
3589  dshape(6,j) = DL0[j];
3590  dshape(8,j) = DL1[j];
3591  dshape(9,j) = DL2[j];
3592  dshape(3,j) = DL3[j];
3593  }
3594  }
3595  else if ((L4 <= 1.0) && (L5 <= 1.0)) // T4
3596  {
3597  for (j = 0; j < 3; j++)
3598  {
3599  dshape(4,j) = - DL5[j];
3600  dshape(5,j) = DL2[j];
3601  dshape(6,j) = - DL4[j];
3602  dshape(8,j) = - DL0[j];
3603  }
3604  }
3605  else if ((L4 >= 1.0) && (L5 <= 1.0)) // T5
3606  {
3607  for (j = 0; j < 3; j++)
3608  {
3609  dshape(4,j) = - DL5[j];
3610  dshape(5,j) = - DL1[j];
3611  dshape(7,j) = DL4[j];
3612  dshape(8,j) = DL3[j];
3613  }
3614  }
3615  else if ((L4 <= 1.0) && (L5 >= 1.0)) // T6
3616  {
3617  for (j = 0; j < 3; j++)
3618  {
3619  dshape(5,j) = - DL3[j];
3620  dshape(6,j) = - DL4[j];
3621  dshape(8,j) = DL1[j];
3622  dshape(9,j) = DL5[j];
3623  }
3624  }
3625  else if ((L4 >= 1.0) && (L5 >= 1.0)) // T7
3626  {
3627  for (j = 0; j < 3; j++)
3628  {
3629  dshape(5,j) = DL0[j];
3630  dshape(7,j) = DL4[j];
3631  dshape(8,j) = - DL2[j];
3632  dshape(9,j) = DL5[j];
3633  }
3634  }
3635 }
3636 
3637 
3638 RefinedBiLinear2DFiniteElement::RefinedBiLinear2DFiniteElement()
3639  : NodalFiniteElement(2, Geometry::SQUARE, 9, 1, FunctionSpace::rQk)
3640 {
3641  Nodes.IntPoint(0).x = 0.0;
3642  Nodes.IntPoint(0).y = 0.0;
3643  Nodes.IntPoint(1).x = 1.0;
3644  Nodes.IntPoint(1).y = 0.0;
3645  Nodes.IntPoint(2).x = 1.0;
3646  Nodes.IntPoint(2).y = 1.0;
3647  Nodes.IntPoint(3).x = 0.0;
3648  Nodes.IntPoint(3).y = 1.0;
3649  Nodes.IntPoint(4).x = 0.5;
3650  Nodes.IntPoint(4).y = 0.0;
3651  Nodes.IntPoint(5).x = 1.0;
3652  Nodes.IntPoint(5).y = 0.5;
3653  Nodes.IntPoint(6).x = 0.5;
3654  Nodes.IntPoint(6).y = 1.0;
3655  Nodes.IntPoint(7).x = 0.0;
3656  Nodes.IntPoint(7).y = 0.5;
3657  Nodes.IntPoint(8).x = 0.5;
3658  Nodes.IntPoint(8).y = 0.5;
3659 }
3660 
3661 void RefinedBiLinear2DFiniteElement::CalcShape(const IntegrationPoint &ip,
3662  Vector &shape) const
3663 {
3664  int i;
3665  double x = ip.x, y = ip.y;
3666  double Lx, Ly;
3667  Lx = 2.0 * ( 1. - x );
3668  Ly = 2.0 * ( 1. - y );
3669 
3670  // The reference square is split in 4 squares as follows:
3671  //
3672  // T0 - 0,4,7,8
3673  // T1 - 1,4,5,8
3674  // T2 - 2,5,6,8
3675  // T3 - 3,6,7,8
3676 
3677  for (i = 0; i < 9; i++)
3678  {
3679  shape(i) = 0.0;
3680  }
3681 
3682  if ((x <= 0.5) && (y <= 0.5)) // T0
3683  {
3684  shape(0) = (Lx - 1.0) * (Ly - 1.0);
3685  shape(4) = (2.0 - Lx) * (Ly - 1.0);
3686  shape(8) = (2.0 - Lx) * (2.0 - Ly);
3687  shape(7) = (Lx - 1.0) * (2.0 - Ly);
3688  }
3689  else if ((x >= 0.5) && (y <= 0.5)) // T1
3690  {
3691  shape(4) = Lx * (Ly - 1.0);
3692  shape(1) = (1.0 - Lx) * (Ly - 1.0);
3693  shape(5) = (1.0 - Lx) * (2.0 - Ly);
3694  shape(8) = Lx * (2.0 - Ly);
3695  }
3696  else if ((x >= 0.5) && (y >= 0.5)) // T2
3697  {
3698  shape(8) = Lx * Ly ;
3699  shape(5) = (1.0 - Lx) * Ly ;
3700  shape(2) = (1.0 - Lx) * (1.0 - Ly);
3701  shape(6) = Lx * (1.0 - Ly);
3702  }
3703  else if ((x <= 0.5) && (y >= 0.5)) // T3
3704  {
3705  shape(7) = (Lx - 1.0) * Ly ;
3706  shape(8) = (2.0 - Lx) * Ly ;
3707  shape(6) = (2.0 - Lx) * (1.0 - Ly);
3708  shape(3) = (Lx - 1.0) * (1.0 - Ly);
3709  }
3710 }
3711 
3712 void RefinedBiLinear2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
3713  DenseMatrix &dshape) const
3714 {
3715  int i,j;
3716  double x = ip.x, y = ip.y;
3717  double Lx, Ly;
3718  Lx = 2.0 * ( 1. - x );
3719  Ly = 2.0 * ( 1. - y );
3720 
3721  for (i = 0; i < 9; i++)
3722  for (j = 0; j < 2; j++)
3723  {
3724  dshape(i,j) = 0.0;
3725  }
3726 
3727  if ((x <= 0.5) && (y <= 0.5)) // T0
3728  {
3729  dshape(0,0) = 2.0 * (1.0 - Ly);
3730  dshape(0,1) = 2.0 * (1.0 - Lx);
3731 
3732  dshape(4,0) = 2.0 * (Ly - 1.0);
3733  dshape(4,1) = -2.0 * (2.0 - Lx);
3734 
3735  dshape(8,0) = 2.0 * (2.0 - Ly);
3736  dshape(8,1) = 2.0 * (2.0 - Lx);
3737 
3738  dshape(7,0) = -2.0 * (2.0 - Ly);
3739  dshape(7,0) = 2.0 * (Lx - 1.0);
3740  }
3741  else if ((x >= 0.5) && (y <= 0.5)) // T1
3742  {
3743  dshape(4,0) = -2.0 * (Ly - 1.0);
3744  dshape(4,1) = -2.0 * Lx;
3745 
3746  dshape(1,0) = 2.0 * (Ly - 1.0);
3747  dshape(1,1) = -2.0 * (1.0 - Lx);
3748 
3749  dshape(5,0) = 2.0 * (2.0 - Ly);
3750  dshape(5,1) = 2.0 * (1.0 - Lx);
3751 
3752  dshape(8,0) = -2.0 * (2.0 - Ly);
3753  dshape(8,1) = 2.0 * Lx;
3754  }
3755  else if ((x >= 0.5) && (y >= 0.5)) // T2
3756  {
3757  dshape(8,0) = -2.0 * Ly;
3758  dshape(8,1) = -2.0 * Lx;
3759 
3760  dshape(5,0) = 2.0 * Ly;
3761  dshape(5,1) = -2.0 * (1.0 - Lx);
3762 
3763  dshape(2,0) = 2.0 * (1.0 - Ly);
3764  dshape(2,1) = 2.0 * (1.0 - Lx);
3765 
3766  dshape(6,0) = -2.0 * (1.0 - Ly);
3767  dshape(6,1) = 2.0 * Lx;
3768  }
3769  else if ((x <= 0.5) && (y >= 0.5)) // T3
3770  {
3771  dshape(7,0) = -2.0 * Ly;
3772  dshape(7,1) = -2.0 * (Lx - 1.0);
3773 
3774  dshape(8,0) = 2.0 * Ly ;
3775  dshape(8,1) = -2.0 * (2.0 - Lx);
3776 
3777  dshape(6,0) = 2.0 * (1.0 - Ly);
3778  dshape(6,1) = 2.0 * (2.0 - Lx);
3779 
3780  dshape(3,0) = -2.0 * (1.0 - Ly);
3781  dshape(3,1) = 2.0 * (Lx - 1.0);
3782  }
3783 }
3784 
3785 RefinedTriLinear3DFiniteElement::RefinedTriLinear3DFiniteElement()
3786  : NodalFiniteElement(3, Geometry::CUBE, 27, 2, FunctionSpace::rQk)
3787 {
3788  double I[27];
3789  double J[27];
3790  double K[27];
3791  // nodes
3792  I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
3793  I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
3794  I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
3795  I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
3796  I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
3797  I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
3798  I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
3799  I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
3800  // edges
3801  I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
3802  I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
3803  I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
3804  I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
3805  I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
3806  I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
3807  I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
3808  I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
3809  I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
3810  I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
3811  I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
3812  I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
3813  // faces
3814  I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
3815  I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
3816  I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
3817  I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
3818  I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
3819  I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
3820  // element
3821  I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
3822 
3823  for (int n = 0; n < 27; n++)
3824  {
3825  Nodes.IntPoint(n).x = I[n];
3826  Nodes.IntPoint(n).y = J[n];
3827  Nodes.IntPoint(n).z = K[n];
3828  }
3829 }
3830 
3831 void RefinedTriLinear3DFiniteElement::CalcShape(const IntegrationPoint &ip,
3832  Vector &shape) const
3833 {
3834  int i, N[8];
3835  double Lx, Ly, Lz;
3836  double x = ip.x, y = ip.y, z = ip.z;
3837 
3838  for (i = 0; i < 27; i++)
3839  {
3840  shape(i) = 0.0;
3841  }
3842 
3843  if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5)) // T0
3844  {
3845  Lx = 1.0 - 2.0 * x;
3846  Ly = 1.0 - 2.0 * y;
3847  Lz = 1.0 - 2.0 * z;
3848 
3849  N[0] = 0;
3850  N[1] = 8;
3851  N[2] = 20;
3852  N[3] = 11;
3853  N[4] = 16;
3854  N[5] = 21;
3855  N[6] = 26;
3856  N[7] = 24;
3857  }
3858  else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5)) // T1
3859  {
3860  Lx = 2.0 - 2.0 * x;
3861  Ly = 1.0 - 2.0 * y;
3862  Lz = 1.0 - 2.0 * z;
3863 
3864  N[0] = 8;
3865  N[1] = 1;
3866  N[2] = 9;
3867  N[3] = 20;
3868  N[4] = 21;
3869  N[5] = 17;
3870  N[6] = 22;
3871  N[7] = 26;
3872  }
3873  else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5)) // T2
3874  {
3875  Lx = 2.0 - 2.0 * x;
3876  Ly = 2.0 - 2.0 * y;
3877  Lz = 1.0 - 2.0 * z;
3878 
3879  N[0] = 20;
3880  N[1] = 9;
3881  N[2] = 2;
3882  N[3] = 10;
3883  N[4] = 26;
3884  N[5] = 22;
3885  N[6] = 18;
3886  N[7] = 23;
3887  }
3888  else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5)) // T3
3889  {
3890  Lx = 1.0 - 2.0 * x;
3891  Ly = 2.0 - 2.0 * y;
3892  Lz = 1.0 - 2.0 * z;
3893 
3894  N[0] = 11;
3895  N[1] = 20;
3896  N[2] = 10;
3897  N[3] = 3;
3898  N[4] = 24;
3899  N[5] = 26;
3900  N[6] = 23;
3901  N[7] = 19;
3902  }
3903  else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5)) // T4
3904  {
3905  Lx = 1.0 - 2.0 * x;
3906  Ly = 1.0 - 2.0 * y;
3907  Lz = 2.0 - 2.0 * z;
3908 
3909  N[0] = 16;
3910  N[1] = 21;
3911  N[2] = 26;
3912  N[3] = 24;
3913  N[4] = 4;
3914  N[5] = 12;
3915  N[6] = 25;
3916  N[7] = 15;
3917  }
3918  else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5)) // T5
3919  {
3920  Lx = 2.0 - 2.0 * x;
3921  Ly = 1.0 - 2.0 * y;
3922  Lz = 2.0 - 2.0 * z;
3923 
3924  N[0] = 21;
3925  N[1] = 17;
3926  N[2] = 22;
3927  N[3] = 26;
3928  N[4] = 12;
3929  N[5] = 5;
3930  N[6] = 13;
3931  N[7] = 25;
3932  }
3933  else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5)) // T6
3934  {
3935  Lx = 2.0 - 2.0 * x;
3936  Ly = 2.0 - 2.0 * y;
3937  Lz = 2.0 - 2.0 * z;
3938 
3939  N[0] = 26;
3940  N[1] = 22;
3941  N[2] = 18;
3942  N[3] = 23;
3943  N[4] = 25;
3944  N[5] = 13;
3945  N[6] = 6;
3946  N[7] = 14;
3947  }
3948  else // T7
3949  {
3950  Lx = 1.0 - 2.0 * x;
3951  Ly = 2.0 - 2.0 * y;
3952  Lz = 2.0 - 2.0 * z;
3953 
3954  N[0] = 24;
3955  N[1] = 26;
3956  N[2] = 23;
3957  N[3] = 19;
3958  N[4] = 15;
3959  N[5] = 25;
3960  N[6] = 14;
3961  N[7] = 7;
3962  }
3963 
3964  shape(N[0]) = Lx * Ly * Lz;
3965  shape(N[1]) = (1 - Lx) * Ly * Lz;
3966  shape(N[2]) = (1 - Lx) * (1 - Ly) * Lz;
3967  shape(N[3]) = Lx * (1 - Ly) * Lz;
3968  shape(N[4]) = Lx * Ly * (1 - Lz);
3969  shape(N[5]) = (1 - Lx) * Ly * (1 - Lz);
3970  shape(N[6]) = (1 - Lx) * (1 - Ly) * (1 - Lz);
3971  shape(N[7]) = Lx * (1 - Ly) * (1 - Lz);
3972 }
3973 
3974 void RefinedTriLinear3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
3975  DenseMatrix &dshape) const
3976 {
3977  int i, j, N[8];
3978  double Lx, Ly, Lz;
3979  double x = ip.x, y = ip.y, z = ip.z;
3980 
3981  for (i = 0; i < 27; i++)
3982  for (j = 0; j < 3; j++)
3983  {
3984  dshape(i,j) = 0.0;
3985  }
3986 
3987  if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5)) // T0
3988  {
3989  Lx = 1.0 - 2.0 * x;
3990  Ly = 1.0 - 2.0 * y;
3991  Lz = 1.0 - 2.0 * z;
3992 
3993  N[0] = 0;
3994  N[1] = 8;
3995  N[2] = 20;
3996  N[3] = 11;
3997  N[4] = 16;
3998  N[5] = 21;
3999  N[6] = 26;
4000  N[7] = 24;
4001  }
4002  else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5)) // T1
4003  {
4004  Lx = 2.0 - 2.0 * x;
4005  Ly = 1.0 - 2.0 * y;
4006  Lz = 1.0 - 2.0 * z;
4007 
4008  N[0] = 8;
4009  N[1] = 1;
4010  N[2] = 9;
4011  N[3] = 20;
4012  N[4] = 21;
4013  N[5] = 17;
4014  N[6] = 22;
4015  N[7] = 26;
4016  }
4017  else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5)) // T2
4018  {
4019  Lx = 2.0 - 2.0 * x;
4020  Ly = 2.0 - 2.0 * y;
4021  Lz = 1.0 - 2.0 * z;
4022 
4023  N[0] = 20;
4024  N[1] = 9;
4025  N[2] = 2;
4026  N[3] = 10;
4027  N[4] = 26;
4028  N[5] = 22;
4029  N[6] = 18;
4030  N[7] = 23;
4031  }
4032  else if ((x >= 0.5) && (y >= 0.5) && (z <= 0.5)) // T3
4033  {
4034  Lx = 1.0 - 2.0 * x;
4035  Ly = 2.0 - 2.0 * y;
4036  Lz = 1.0 - 2.0 * z;
4037 
4038  N[0] = 11;
4039  N[1] = 20;
4040  N[2] = 10;
4041  N[3] = 3;
4042  N[4] = 24;
4043  N[5] = 26;
4044  N[6] = 23;
4045  N[7] = 19;
4046  }
4047  else if ((x <= 0.5) && (y <= 0.5) && (z >= 0.5)) // T4
4048  {
4049  Lx = 1.0 - 2.0 * x;
4050  Ly = 1.0 - 2.0 * y;
4051  Lz = 2.0 - 2.0 * z;
4052 
4053  N[0] = 16;
4054  N[1] = 21;
4055  N[2] = 26;
4056  N[3] = 24;
4057  N[4] = 4;
4058  N[5] = 12;
4059  N[6] = 25;
4060  N[7] = 15;
4061  }
4062  else if ((x >= 0.5) && (y <= 0.5) && (z >= 0.5)) // T5
4063  {
4064  Lx = 2.0 - 2.0 * x;
4065  Ly = 1.0 - 2.0 * y;
4066  Lz = 2.0 - 2.0 * z;
4067 
4068  N[0] = 21;
4069  N[1] = 17;
4070  N[2] = 22;
4071  N[3] = 26;
4072  N[4] = 12;
4073  N[5] = 5;
4074  N[6] = 13;
4075  N[7] = 25;
4076  }
4077  else if ((x <= 0.5) && (y >= 0.5) && (z >= 0.5)) // T6
4078  {
4079  Lx = 2.0 - 2.0 * x;
4080  Ly = 2.0 - 2.0 * y;
4081  Lz = 2.0 - 2.0 * z;
4082 
4083  N[0] = 26;
4084  N[1] = 22;
4085  N[2] = 18;
4086  N[3] = 23;
4087  N[4] = 25;
4088  N[5] = 13;
4089  N[6] = 6;
4090  N[7] = 14;
4091  }
4092  else // T7
4093  {
4094  Lx = 1.0 - 2.0 * x;
4095  Ly = 2.0 - 2.0 * y;
4096  Lz = 2.0 - 2.0 * z;
4097 
4098  N[0] = 24;
4099  N[1] = 26;
4100  N[2] = 23;
4101  N[3] = 19;
4102  N[4] = 15;
4103  N[5] = 25;
4104  N[6] = 14;
4105  N[7] = 7;
4106  }
4107 
4108  dshape(N[0],0) = -2.0 * Ly * Lz ;
4109  dshape(N[0],1) = -2.0 * Lx * Lz ;
4110  dshape(N[0],2) = -2.0 * Lx * Ly ;
4111 
4112  dshape(N[1],0) = 2.0 * Ly * Lz ;
4113  dshape(N[1],1) = -2.0 * (1 - Lx) * Lz ;
4114  dshape(N[1],2) = -2.0 * (1 - Lx) * Ly ;
4115 
4116  dshape(N[2],0) = 2.0 * (1 - Ly) * Lz ;
4117  dshape(N[2],1) = 2.0 * (1 - Lx) * Lz ;
4118  dshape(N[2],2) = -2.0 * (1 - Lx) * (1 - Ly);
4119 
4120  dshape(N[3],0) = -2.0 * (1 - Ly) * Lz ;
4121  dshape(N[3],1) = 2.0 * Lx * Lz ;
4122  dshape(N[3],2) = -2.0 * Lx * (1 - Ly);
4123 
4124  dshape(N[4],0) = -2.0 * Ly * (1 - Lz);
4125  dshape(N[4],1) = -2.0 * Lx * (1 - Lz);
4126  dshape(N[4],2) = 2.0 * Lx * Ly ;
4127 
4128  dshape(N[5],0) = 2.0 * Ly * (1 - Lz);
4129  dshape(N[5],1) = -2.0 * (1 - Lx) * (1 - Lz);
4130  dshape(N[5],2) = 2.0 * (1 - Lx) * Ly ;
4131 
4132  dshape(N[6],0) = 2.0 * (1 - Ly) * (1 - Lz);
4133  dshape(N[6],1) = 2.0 * (1 - Lx) * (1 - Lz);
4134  dshape(N[6],2) = 2.0 * (1 - Lx) * (1 - Ly);
4135 
4136  dshape(N[7],0) = -2.0 * (1 - Ly) * (1 - Lz);
4137  dshape(N[7],1) = 2.0 * Lx * (1 - Lz);
4138  dshape(N[7],2) = 2.0 * Lx * (1 - Ly);
4139 }
4140 
4141 
4142 Nedelec1HexFiniteElement::Nedelec1HexFiniteElement()
4143  : VectorFiniteElement(3, Geometry::CUBE, 12, 1, H_CURL,
4144  FunctionSpace::Qk)
4145 {
4146  // not real nodes ...
4147  Nodes.IntPoint(0).x = 0.5;
4148  Nodes.IntPoint(0).y = 0.0;
4149  Nodes.IntPoint(0).z = 0.0;
4150 
4151  Nodes.IntPoint(1).x = 1.0;
4152  Nodes.IntPoint(1).y = 0.5;
4153  Nodes.IntPoint(1).z = 0.0;
4154 
4155  Nodes.IntPoint(2).x = 0.5;
4156  Nodes.IntPoint(2).y = 1.0;
4157  Nodes.IntPoint(2).z = 0.0;
4158 
4159  Nodes.IntPoint(3).x = 0.0;
4160  Nodes.IntPoint(3).y = 0.5;
4161  Nodes.IntPoint(3).z = 0.0;
4162 
4163  Nodes.IntPoint(4).x = 0.5;
4164  Nodes.IntPoint(4).y = 0.0;
4165  Nodes.IntPoint(4).z = 1.0;
4166 
4167  Nodes.IntPoint(5).x = 1.0;
4168  Nodes.IntPoint(5).y = 0.5;
4169  Nodes.IntPoint(5).z = 1.0;
4170 
4171  Nodes.IntPoint(6).x = 0.5;
4172  Nodes.IntPoint(6).y = 1.0;
4173  Nodes.IntPoint(6).z = 1.0;
4174 
4175  Nodes.IntPoint(7).x = 0.0;
4176  Nodes.IntPoint(7).y = 0.5;
4177  Nodes.IntPoint(7).z = 1.0;
4178 
4179  Nodes.IntPoint(8).x = 0.0;
4180  Nodes.IntPoint(8).y = 0.0;
4181  Nodes.IntPoint(8).z = 0.5;
4182 
4183  Nodes.IntPoint(9).x = 1.0;
4184  Nodes.IntPoint(9).y = 0.0;
4185  Nodes.IntPoint(9).z = 0.5;
4186 
4187  Nodes.IntPoint(10).x= 1.0;
4188  Nodes.IntPoint(10).y= 1.0;
4189  Nodes.IntPoint(10).z= 0.5;
4190 
4191  Nodes.IntPoint(11).x= 0.0;
4192  Nodes.IntPoint(11).y= 1.0;
4193  Nodes.IntPoint(11).z= 0.5;
4194 }
4195 
4196 void Nedelec1HexFiniteElement::CalcVShape(const IntegrationPoint &ip,
4197  DenseMatrix &shape) const
4198 {
4199  double x = ip.x, y = ip.y, z = ip.z;
4200 
4201  shape(0,0) = (1. - y) * (1. - z);
4202  shape(0,1) = 0.;
4203  shape(0,2) = 0.;
4204 
4205  shape(2,0) = y * (1. - z);
4206  shape(2,1) = 0.;
4207  shape(2,2) = 0.;
4208 
4209  shape(4,0) = z * (1. - y);
4210  shape(4,1) = 0.;
4211  shape(4,2) = 0.;
4212 
4213  shape(6,0) = y * z;
4214  shape(6,1) = 0.;
4215  shape(6,2) = 0.;
4216 
4217  shape(1,0) = 0.;
4218  shape(1,1) = x * (1. - z);
4219  shape(1,2) = 0.;
4220 
4221  shape(3,0) = 0.;
4222  shape(3,1) = (1. - x) * (1. - z);
4223  shape(3,2) = 0.;
4224 
4225  shape(5,0) = 0.;
4226  shape(5,1) = x * z;
4227  shape(5,2) = 0.;
4228 
4229  shape(7,0) = 0.;
4230  shape(7,1) = (1. - x) * z;
4231  shape(7,2) = 0.;
4232 
4233  shape(8,0) = 0.;
4234  shape(8,1) = 0.;
4235  shape(8,2) = (1. - x) * (1. - y);
4236 
4237  shape(9,0) = 0.;
4238  shape(9,1) = 0.;
4239  shape(9,2) = x * (1. - y);
4240 
4241  shape(10,0) = 0.;
4242  shape(10,1) = 0.;
4243  shape(10,2) = x * y;
4244 
4245  shape(11,0) = 0.;
4246  shape(11,1) = 0.;
4247  shape(11,2) = y * (1. - x);
4248 
4249 }
4250 
4251 void Nedelec1HexFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
4252  DenseMatrix &curl_shape)
4253 const
4254 {
4255  double x = ip.x, y = ip.y, z = ip.z;
4256 
4257  curl_shape(0,0) = 0.;
4258  curl_shape(0,1) = y - 1.;
4259  curl_shape(0,2) = 1. - z;
4260 
4261  curl_shape(2,0) = 0.;
4262  curl_shape(2,1) = -y;
4263  curl_shape(2,2) = z - 1.;
4264 
4265  curl_shape(4,0) = 0;
4266  curl_shape(4,1) = 1. - y;
4267  curl_shape(4,2) = z;
4268 
4269  curl_shape(6,0) = 0.;
4270  curl_shape(6,1) = y;
4271  curl_shape(6,2) = -z;
4272 
4273  curl_shape(1,0) = x;
4274  curl_shape(1,1) = 0.;
4275  curl_shape(1,2) = 1. - z;
4276 
4277  curl_shape(3,0) = 1. - x;
4278  curl_shape(3,1) = 0.;
4279  curl_shape(3,2) = z - 1.;
4280 
4281  curl_shape(5,0) = -x;
4282  curl_shape(5,1) = 0.;
4283  curl_shape(5,2) = z;
4284 
4285  curl_shape(7,0) = x - 1.;
4286  curl_shape(7,1) = 0.;
4287  curl_shape(7,2) = -z;
4288 
4289  curl_shape(8,0) = x - 1.;
4290  curl_shape(8,1) = 1. - y;
4291  curl_shape(8,2) = 0.;
4292 
4293  curl_shape(9,0) = -x;
4294  curl_shape(9,1) = y - 1.;
4295  curl_shape(9,2) = 0;
4296 
4297  curl_shape(10,0) = x;
4298  curl_shape(10,1) = -y;
4299  curl_shape(10,2) = 0.;
4300 
4301  curl_shape(11,0) = 1. - x;
4302  curl_shape(11,1) = y;
4303  curl_shape(11,2) = 0.;
4304 }
4305 
4306 const double Nedelec1HexFiniteElement::tk[12][3] =
4307 {
4308  {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4309  {1,0,0}, {0,1,0}, {1,0,0}, {0,1,0},
4310  {0,0,1}, {0,0,1}, {0,0,1}, {0,0,1}
4311 };
4312 
4313 void Nedelec1HexFiniteElement::GetLocalInterpolation (
4314  ElementTransformation &Trans, DenseMatrix &I) const
4315 {
4316  int k, j;
4317 #ifdef MFEM_THREAD_SAFE
4318  DenseMatrix vshape(dof, dim);
4319 #endif
4320 
4321 #ifdef MFEM_DEBUG
4322  for (k = 0; k < dof; k++)
4323  {
4324  CalcVShape (Nodes.IntPoint(k), vshape);
4325  for (j = 0; j < dof; j++)
4326  {
4327  double d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
4328  vshape(j,2)*tk[k][2] );
4329  if (j == k) { d -= 1.0; }
4330  if (fabs(d) > 1.0e-12)
4331  {
4332  mfem::err << "Nedelec1HexFiniteElement::GetLocalInterpolation (...)\n"
4333  " k = " << k << ", j = " << j << ", d = " << d << endl;
4334  mfem_error();
4335  }
4336  }
4337  }
4338 #endif
4339 
4340  IntegrationPoint ip;
4341  ip.x = ip.y = ip.z = 0.0;
4342  Trans.SetIntPoint (&ip);
4343  // Trans must be linear (more to have embedding?)
4344  const DenseMatrix &J = Trans.Jacobian();
4345  double vk[3];
4346  Vector xk (vk, 3);
4347 
4348  for (k = 0; k < dof; k++)
4349  {
4350  Trans.Transform (Nodes.IntPoint (k), xk);
4351  ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
4352  CalcVShape (ip, vshape);
4353  // vk = J tk
4354  vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4355  vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4356  vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4357  for (j = 0; j < dof; j++)
4358  if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
4359  vshape(j,2)*vk[2])) < 1.0e-12)
4360  {
4361  I(k,j) = 0.0;
4362  }
4363  }
4364 }
4365 
4366 void Nedelec1HexFiniteElement::Project (
4367  VectorCoefficient &vc, ElementTransformation &Trans,
4368  Vector &dofs) const
4369 {
4370  double vk[3];
4371  Vector xk (vk, 3);
4372 
4373  for (int k = 0; k < dof; k++)
4374  {
4375  Trans.SetIntPoint (&Nodes.IntPoint (k));
4376  const DenseMatrix &J = Trans.Jacobian();
4377 
4378  vc.Eval (xk, Trans, Nodes.IntPoint (k));
4379  // xk^t J tk
4380  dofs(k) =
4381  vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4382  vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4383  vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4384  }
4385 }
4386 
4387 void Nedelec1HexFiniteElement::ProjectGrad(const FiniteElement &fe,
4388  ElementTransformation &Trans,
4389  DenseMatrix &grad) const
4390 {
4391  DenseMatrix dshape(fe.GetDof(), 3);
4392  Vector grad_k(fe.GetDof());
4393 
4394  grad.SetSize(dof, fe.GetDof());
4395  for (int k = 0; k < dof; k++)
4396  {
4397  fe.CalcDShape(Nodes.IntPoint(k), dshape);
4398  dshape.Mult(tk[k], grad_k);
4399  for (int j = 0; j < grad_k.Size(); j++)
4400  {
4401  grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4402  }
4403  }
4404 }
4405 
4406 
4407 Nedelec1TetFiniteElement::Nedelec1TetFiniteElement()
4408  : VectorFiniteElement(3, Geometry::TETRAHEDRON, 6, 1, H_CURL)
4409 {
4410  // not real nodes ...
4411  Nodes.IntPoint(0).x = 0.5;
4412  Nodes.IntPoint(0).y = 0.0;
4413  Nodes.IntPoint(0).z = 0.0;
4414 
4415  Nodes.IntPoint(1).x = 0.0;
4416  Nodes.IntPoint(1).y = 0.5;
4417  Nodes.IntPoint(1).z = 0.0;
4418 
4419  Nodes.IntPoint(2).x = 0.0;
4420  Nodes.IntPoint(2).y = 0.0;
4421  Nodes.IntPoint(2).z = 0.5;
4422 
4423  Nodes.IntPoint(3).x = 0.5;
4424  Nodes.IntPoint(3).y = 0.5;
4425  Nodes.IntPoint(3).z = 0.0;
4426 
4427  Nodes.IntPoint(4).x = 0.5;
4428  Nodes.IntPoint(4).y = 0.0;
4429  Nodes.IntPoint(4).z = 0.5;
4430 
4431  Nodes.IntPoint(5).x = 0.0;
4432  Nodes.IntPoint(5).y = 0.5;
4433  Nodes.IntPoint(5).z = 0.5;
4434 }
4435 
4436 void Nedelec1TetFiniteElement::CalcVShape(const IntegrationPoint &ip,
4437  DenseMatrix &shape) const
4438 {
4439  double x = ip.x, y = ip.y, z = ip.z;
4440 
4441  shape(0,0) = 1. - y - z;
4442  shape(0,1) = x;
4443  shape(0,2) = x;
4444 
4445  shape(1,0) = y;
4446  shape(1,1) = 1. - x - z;
4447  shape(1,2) = y;
4448 
4449  shape(2,0) = z;
4450  shape(2,1) = z;
4451  shape(2,2) = 1. - x - y;
4452 
4453  shape(3,0) = -y;
4454  shape(3,1) = x;
4455  shape(3,2) = 0.;
4456 
4457  shape(4,0) = -z;
4458  shape(4,1) = 0.;
4459  shape(4,2) = x;
4460 
4461  shape(5,0) = 0.;
4462  shape(5,1) = -z;
4463  shape(5,2) = y;
4464 }
4465 
4466 void Nedelec1TetFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
4467  DenseMatrix &curl_shape)
4468 const
4469 {
4470  curl_shape(0,0) = 0.;
4471  curl_shape(0,1) = -2.;
4472  curl_shape(0,2) = 2.;
4473 
4474  curl_shape(1,0) = 2.;
4475  curl_shape(1,1) = 0.;
4476  curl_shape(1,2) = -2.;
4477 
4478  curl_shape(2,0) = -2.;
4479  curl_shape(2,1) = 2.;
4480  curl_shape(2,2) = 0.;
4481 
4482  curl_shape(3,0) = 0.;
4483  curl_shape(3,1) = 0.;
4484  curl_shape(3,2) = 2.;
4485 
4486  curl_shape(4,0) = 0.;
4487  curl_shape(4,1) = -2.;
4488  curl_shape(4,2) = 0.;
4489 
4490  curl_shape(5,0) = 2.;
4491  curl_shape(5,1) = 0.;
4492  curl_shape(5,2) = 0.;
4493 }
4494 
4495 const double Nedelec1TetFiniteElement::tk[6][3] =
4496 {{1,0,0}, {0,1,0}, {0,0,1}, {-1,1,0}, {-1,0,1}, {0,-1,1}};
4497 
4498 void Nedelec1TetFiniteElement::GetLocalInterpolation (
4499  ElementTransformation &Trans, DenseMatrix &I) const
4500 {
4501  int k, j;
4502 #ifdef MFEM_THREAD_SAFE
4503  DenseMatrix vshape(dof, dim);
4504 #endif
4505 
4506 #ifdef MFEM_DEBUG
4507  for (k = 0; k < dof; k++)
4508  {
4509  CalcVShape (Nodes.IntPoint(k), vshape);
4510  for (j = 0; j < dof; j++)
4511  {
4512  double d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
4513  vshape(j,2)*tk[k][2] );
4514  if (j == k) { d -= 1.0; }
4515  if (fabs(d) > 1.0e-12)
4516  {
4517  mfem::err << "Nedelec1TetFiniteElement::GetLocalInterpolation (...)\n"
4518  " k = " << k << ", j = " << j << ", d = " << d << endl;
4519  mfem_error();
4520  }
4521  }
4522  }
4523 #endif
4524 
4525  IntegrationPoint ip;
4526  ip.x = ip.y = ip.z = 0.0;
4527  Trans.SetIntPoint (&ip);
4528  // Trans must be linear
4529  const DenseMatrix &J = Trans.Jacobian();
4530  double vk[3];
4531  Vector xk (vk, 3);
4532 
4533  for (k = 0; k < dof; k++)
4534  {
4535  Trans.Transform (Nodes.IntPoint (k), xk);
4536  ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
4537  CalcVShape (ip, vshape);
4538  // vk = J tk
4539  vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4540  vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4541  vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4542  for (j = 0; j < dof; j++)
4543  if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
4544  vshape(j,2)*vk[2])) < 1.0e-12)
4545  {
4546  I(k,j) = 0.0;
4547  }
4548  }
4549 }
4550 
4551 void Nedelec1TetFiniteElement::Project (
4552  VectorCoefficient &vc, ElementTransformation &Trans,
4553  Vector &dofs) const
4554 {
4555  double vk[3];
4556  Vector xk (vk, 3);
4557 
4558  for (int k = 0; k < dof; k++)
4559  {
4560  Trans.SetIntPoint (&Nodes.IntPoint (k));
4561  const DenseMatrix &J = Trans.Jacobian();
4562 
4563  vc.Eval (xk, Trans, Nodes.IntPoint (k));
4564  // xk^t J tk
4565  dofs(k) =
4566  vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4567  vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4568  vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4569  }
4570 }
4571 
4572 void Nedelec1TetFiniteElement::ProjectGrad(const FiniteElement &fe,
4573  ElementTransformation &Trans,
4574  DenseMatrix &grad) const
4575 {
4576  DenseMatrix dshape(fe.GetDof(), 3);
4577  Vector grad_k(fe.GetDof());
4578 
4579  grad.SetSize(dof, fe.GetDof());
4580  for (int k = 0; k < dof; k++)
4581  {
4582  fe.CalcDShape(Nodes.IntPoint(k), dshape);
4583  dshape.Mult(tk[k], grad_k);
4584  for (int j = 0; j < grad_k.Size(); j++)
4585  {
4586  grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4587  }
4588  }
4589 }
4590 
4591 
4592 Nedelec1WdgFiniteElement::Nedelec1WdgFiniteElement()
4593  : VectorFiniteElement(3, Geometry::PRISM, 9, 1, H_CURL, FunctionSpace::Qk)
4594 {
4595  // not real nodes ...
4596  Nodes.IntPoint(0).x = 0.5;
4597  Nodes.IntPoint(0).y = 0.0;
4598  Nodes.IntPoint(0).z = 0.0;
4599 
4600  Nodes.IntPoint(1).x = 0.5;
4601  Nodes.IntPoint(1).y = 0.5;
4602  Nodes.IntPoint(1).z = 0.0;
4603 
4604  Nodes.IntPoint(2).x = 0.0;
4605  Nodes.IntPoint(2).y = 0.5;
4606  Nodes.IntPoint(2).z = 0.0;
4607 
4608  Nodes.IntPoint(3).x = 0.5;
4609  Nodes.IntPoint(3).y = 0.0;
4610  Nodes.IntPoint(3).z = 1.0;
4611 
4612  Nodes.IntPoint(4).x = 0.5;
4613  Nodes.IntPoint(4).y = 0.5;
4614  Nodes.IntPoint(4).z = 1.0;
4615 
4616  Nodes.IntPoint(5).x = 0.0;
4617  Nodes.IntPoint(5).y = 0.5;
4618  Nodes.IntPoint(5).z = 1.0;
4619 
4620  Nodes.IntPoint(6).x = 0.0;
4621  Nodes.IntPoint(6).y = 0.0;
4622  Nodes.IntPoint(6).z = 0.5;
4623 
4624  Nodes.IntPoint(7).x = 1.0;
4625  Nodes.IntPoint(7).y = 0.0;
4626  Nodes.IntPoint(7).z = 0.5;
4627 
4628  Nodes.IntPoint(8).x = 0.0;
4629  Nodes.IntPoint(8).y = 1.0;
4630  Nodes.IntPoint(8).z = 0.5;
4631 }
4632 
4633 void Nedelec1WdgFiniteElement::CalcVShape(const IntegrationPoint &ip,
4634  DenseMatrix &shape) const
4635 {
4636  double x = ip.x, y = ip.y, z = ip.z;
4637 
4638  shape(0,0) = (1. - y) * (1. - z);
4639  shape(0,1) = x * (1. - z);
4640  shape(0,2) = 0.;
4641 
4642  shape(1,0) = - y * (1. - z);
4643  shape(1,1) = x * (1. - z);
4644  shape(1,2) = 0.;
4645 
4646  shape(2,0) = - y * (1. - z);
4647  shape(2,1) = - (1. - x) * (1. - z);
4648  shape(2,2) = 0.;
4649 
4650  shape(3,0) = (1. - y) * z;
4651  shape(3,1) = x * z;
4652  shape(3,2) = 0.;
4653 
4654  shape(4,0) = - y * z;
4655  shape(4,1) = x * z;
4656  shape(4,2) = 0.;
4657 
4658  shape(5,0) = - y * z;
4659  shape(5,1) = - (1. - x) * z;
4660  shape(5,2) = 0.;
4661 
4662  shape(6,0) = 0.;
4663  shape(6,1) = 0.;
4664  shape(6,2) = 1. - x - y;
4665 
4666  shape(7,0) = 0.;
4667  shape(7,1) = 0.;
4668  shape(7,2) = x;
4669 
4670  shape(8,0) = 0.;
4671  shape(8,1) = 0.;
4672  shape(8,2) = y;
4673 }
4674 
4675 void Nedelec1WdgFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
4676  DenseMatrix &curl_shape)
4677 const
4678 {
4679  double x = ip.x, y = ip.y, z2 = 2. * ip.z;
4680 
4681  curl_shape(0,0) = x;
4682  curl_shape(0,1) = - 1. + y;
4683  curl_shape(0,2) = 2. - z2;
4684 
4685  curl_shape(1,0) = x;
4686  curl_shape(1,1) = y;
4687  curl_shape(1,2) = 2. - z2;
4688 
4689  curl_shape(2,0) = - 1. + x;
4690  curl_shape(2,1) = y;
4691  curl_shape(2,2) = 2. - z2;
4692 
4693  curl_shape(3,0) = - x;
4694  curl_shape(3,1) = 1. - y;
4695  curl_shape(3,2) = z2;
4696 
4697  curl_shape(4,0) = - x;
4698  curl_shape(4,1) = - y;
4699  curl_shape(4,2) = z2;
4700 
4701  curl_shape(5,0) = 1. - x;
4702  curl_shape(5,1) = - y;
4703  curl_shape(5,2) = z2;
4704 
4705  curl_shape(6,0) = - 1.;
4706  curl_shape(6,1) = 1.;
4707  curl_shape(6,2) = 0.;
4708 
4709  curl_shape(7,0) = 0.;
4710  curl_shape(7,1) = - 1.;
4711  curl_shape(7,2) = 0.;
4712 
4713  curl_shape(8,0) = 1.;
4714  curl_shape(8,1) = 0.;
4715  curl_shape(8,2) = 0.;
4716 }
4717 
4718 const double Nedelec1WdgFiniteElement::tk[9][3] =
4719 {
4720  {1,0,0}, {-1,1,0}, {0,-1,0}, {1,0,0}, {-1,1,0}, {0,-1,0},
4721  {0,0,1}, {0,0,1}, {0,0,1}
4722 };
4723 
4724 void Nedelec1WdgFiniteElement::GetLocalInterpolation (
4725  ElementTransformation &Trans, DenseMatrix &I) const
4726 {
4727  int k, j;
4728 #ifdef MFEM_THREAD_SAFE
4729  DenseMatrix vshape(dof, dim);
4730 #endif
4731 
4732 #ifdef MFEM_DEBUG
4733  for (k = 0; k < dof; k++)
4734  {
4735  CalcVShape (Nodes.IntPoint(k), vshape);
4736  for (j = 0; j < dof; j++)
4737  {
4738  double d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
4739  vshape(j,2)*tk[k][2] );
4740  if (j == k) { d -= 1.0; }
4741  if (fabs(d) > 1.0e-12)
4742  {
4743  mfem::err << "Nedelec1WdgFiniteElement::GetLocalInterpolation (...)\n"
4744  " k = " << k << ", j = " << j << ", d = " << d << endl;
4745  mfem_error();
4746  }
4747  }
4748  }
4749 #endif
4750 
4751  IntegrationPoint ip;
4752  ip.x = ip.y = ip.z = 0.0;
4753  Trans.SetIntPoint (&ip);
4754  // Trans must be linear
4755  const DenseMatrix &J = Trans.Jacobian();
4756  double vk[3];
4757  Vector xk (vk, 3);
4758 
4759  for (k = 0; k < dof; k++)
4760  {
4761  Trans.Transform (Nodes.IntPoint (k), xk);
4762  ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
4763  CalcVShape (ip, vshape);
4764  // vk = J tk
4765  vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
4766  vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
4767  vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
4768  for (j = 0; j < dof; j++)
4769  if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
4770  vshape(j,2)*vk[2])) < 1.0e-12)
4771  {
4772  I(k,j) = 0.0;
4773  }
4774  }
4775 }
4776 
4777 void Nedelec1WdgFiniteElement::Project (
4778  VectorCoefficient &vc, ElementTransformation &Trans,
4779  Vector &dofs) const
4780 {
4781  double vk[3];
4782  Vector xk (vk, 3);
4783 
4784  for (int k = 0; k < dof; k++)
4785  {
4786  Trans.SetIntPoint (&Nodes.IntPoint (k));
4787  const DenseMatrix &J = Trans.Jacobian();
4788 
4789  vc.Eval (xk, Trans, Nodes.IntPoint (k));
4790  // xk^t J tk
4791  dofs(k) =
4792  vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
4793  vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
4794  vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
4795  }
4796 }
4797 
4798 void Nedelec1WdgFiniteElement::ProjectGrad(const FiniteElement &fe,
4799  ElementTransformation &Trans,
4800  DenseMatrix &grad) const
4801 {
4802  DenseMatrix dshape(fe.GetDof(), 3);
4803  Vector grad_k(fe.GetDof());
4804 
4805  grad.SetSize(dof, fe.GetDof());
4806  for (int k = 0; k < dof; k++)
4807  {
4808  fe.CalcDShape(Nodes.IntPoint(k), dshape);
4809  dshape.Mult(tk[k], grad_k);
4810  for (int j = 0; j < grad_k.Size(); j++)
4811  {
4812  grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
4813  }
4814  }
4815 }
4816 
4817 
4818 Nedelec1PyrFiniteElement::Nedelec1PyrFiniteElement()
4819  : VectorFiniteElement(3, Geometry::PYRAMID, 8, 1, H_CURL)
4820 {
4821  // not real nodes ...
4822  Nodes.IntPoint(0).x = 0.5;
4823  Nodes.IntPoint(0).y = 0.0;
4824  Nodes.IntPoint(0).z = 0.0;
4825 
4826  Nodes.IntPoint(1).x = 1.0;
4827  Nodes.IntPoint(1).y = 0.5;
4828  Nodes.IntPoint(1).z = 0.0;
4829 
4830  Nodes.IntPoint(2).x = 0.5;
4831  Nodes.IntPoint(2).y = 1.0;
4832  Nodes.IntPoint(2).z = 0.0;
4833 
4834  Nodes.IntPoint(3).x = 0.0;
4835  Nodes.IntPoint(3).y = 0.5;
4836  Nodes.IntPoint(3).z = 0.0;
4837 
4838  Nodes.IntPoint(4).x = 0.0;
4839  Nodes.IntPoint(4).y = 0.0;
4840  Nodes.IntPoint(4).z = 0.5;
4841 
4842  Nodes.IntPoint(5).x = 0.5;
4843  Nodes.IntPoint(5).y = 0.0;
4844  Nodes.IntPoint(5).z = 0.5;
4845 
4846  Nodes.IntPoint(6).x = 0.5;
4847  Nodes.IntPoint(6).y = 0.5;
4848  Nodes.IntPoint(6).z = 0.5;
4849 
4850  Nodes.IntPoint(7).x = 0.0;
4851  Nodes.IntPoint(7).y = 0.5;
4852  Nodes.IntPoint(7).z = 0.5;
4853 }
4854 
4855 void Nedelec1PyrFiniteElement::CalcVShape(const IntegrationPoint &ip,
4856  DenseMatrix &shape) const
4857 {
4858  double x = ip.x, y = ip.y, z = ip.z, z2 = 2. * ip.z;
4859  double ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4860 
4861  double tol = 1e-6;
4862 
4863  if (oz <= tol)
4864  {
4865  // We must return the limit of the basis functions as z->1. In order to
4866  // remain inside the pyramid in this limit the x and y coordinates must
4867  // be approaching 0. The resulting limiting basis function values are:
4868  shape(0,0) = 0.;
4869  shape(0,1) = 0.;
4870  shape(0,2) = 0.;
4871 
4872  shape(1,0) = 0.;
4873  shape(1,1) = 0.;
4874  shape(1,2) = 0.;
4875 
4876  shape(2,0) = 0.;
4877  shape(2,1) = 0.;
4878  shape(2,2) = 0.;
4879 
4880  shape(3,0) = 0.;
4881  shape(3,1) = 0.;
4882  shape(3,2) = 0.;
4883 
4884  shape(4,0) = 1.;
4885  shape(4,1) = 1.;
4886  shape(4,2) = 1.;
4887 
4888  shape(5,0) = - 1.;
4889  shape(5,1) = 0.;
4890  shape(5,2) = 0.;
4891 
4892  shape(6,0) = 0.;
4893  shape(6,1) = 0.;
4894  shape(6,2) = 0.;
4895 
4896  shape(7,0) = 0.;
4897  shape(7,1) = - 1.;
4898  shape(7,2) = 0.;
4899 
4900  return;
4901  }
4902 
4903  double ozi = 1.0 / oz;
4904 
4905  shape(0,0) = oy;
4906  shape(0,1) = 0.;
4907  shape(0,2) = x * oy * ozi;
4908 
4909  shape(1,0) = 0.;
4910  shape(1,1) = x;
4911  shape(1,2) = x * y * ozi;
4912 
4913  shape(2,0) = y;
4914  shape(2,1) = 0.;
4915  shape(2,2) = x * y * ozi;
4916 
4917  shape(3,0) = 0.;
4918  shape(3,1) = ox;
4919  shape(3,2) = ox * y * ozi;
4920 
4921  shape(4,0) = oy * z * ozi;
4922  shape(4,1) = ox * z * ozi;
4923  shape(4,2) = 1. - x - y + x * y * (1. - z2) * ozi * ozi;
4924 
4925  shape(5,0) = - oy * z * ozi;
4926  shape(5,1) = x * z * ozi;
4927  shape(5,2) = x * (1. - y * (1. - z2) * ozi * ozi);
4928 
4929  shape(6,0) = - y * z * ozi;
4930  shape(6,1) = - x * z * ozi;
4931  shape(6,2) = x * y * (1. - z2) * ozi * ozi;
4932 
4933  shape(7,0) = y * z * ozi;
4934  shape(7,1) = - ox * z * ozi;
4935  shape(7,2) = y * (1. - x * (1. - z2) * ozi * ozi);
4936 }
4937 
4938 void Nedelec1PyrFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
4939  DenseMatrix &curl_shape)
4940 const
4941 {
4942  double x = ip.x, y = ip.y, z = ip.z, z2 = 2. * z;
4943  double ox = 1. - x - z, oy = 1. - y - z, oz = 1. - z;
4944 
4945  double tol = 1e-6;
4946 
4947  if (oz <= tol)
4948  {
4949  // We must return the limit of the basis function derivatives as z->1.
4950  // In order to remain inside the pyramid in this limit the x and y
4951  // coordinates must be approaching 0. The resulting limiting basis
4952  // function values are:
4953  curl_shape(0,0) = 0.;
4954  curl_shape(0,1) = - 2.;
4955  curl_shape(0,2) = 1.;
4956 
4957  curl_shape(1,0) = 0.;
4958  curl_shape(1,1) = 0.;
4959  curl_shape(1,2) = 1.;
4960 
4961  curl_shape(2,0) = 0.;
4962  curl_shape(2,1) = 0.;
4963  curl_shape(2,2) = - 1.;
4964 
4965  curl_shape(3,0) = 2.;
4966  curl_shape(3,1) = 0.;
4967  curl_shape(3,2) = - 1.;
4968 
4969  curl_shape(4,0) = - 2.;
4970  curl_shape(4,1) = 2.;
4971  curl_shape(4,2) = 0.;
4972 
4973  curl_shape(5,0) = 0.;
4974  curl_shape(5,1) = - 2.;
4975  curl_shape(5,2) = 0.;
4976 
4977  curl_shape(6,0) = 0.;
4978  curl_shape(6,1) = 0.;
4979  curl_shape(6,2) = 0.;
4980 
4981  curl_shape(7,0) = 2.;
4982  curl_shape(7,1) = 0.;
4983  curl_shape(7,2) = 0.;
4984 
4985  return;
4986  }
4987 
4988  double ozi = 1. / oz;
4989 
4990  curl_shape(0,0) = - x * ozi;
4991  curl_shape(0,1) = - 2. + y * ozi;
4992  curl_shape(0,2) = 1.;
4993 
4994  curl_shape(1,0) = x * ozi;
4995  curl_shape(1,1) = - y * ozi;
4996  curl_shape(1,2) = 1.;
4997 
4998  curl_shape(2,0) = x * ozi;
4999  curl_shape(2,1) = - y * ozi;
5000  curl_shape(2,2) = - 1.;
5001 
5002  curl_shape(3,0) = (2. - x - z2) * ozi;
5003  curl_shape(3,1) = y * ozi;
5004  curl_shape(3,2) = - 1.;
5005 
5006  curl_shape(4,0) = - 2. * ox * ozi;
5007  curl_shape(4,1) = 2. * oy * ozi;
5008  curl_shape(4,2) = 0.;
5009 
5010  curl_shape(5,0) = - 2. * x * ozi;
5011  curl_shape(5,1) = - 2. * oy * ozi;
5012  curl_shape(5,2) = 0.;
5013 
5014  curl_shape(6,0) = 2. * x * ozi;
5015  curl_shape(6,1) = - 2. * y * ozi;
5016  curl_shape(6,2) = 0.;
5017 
5018  curl_shape(7,0) = 2. * ox * ozi;
5019  curl_shape(7,1) = 2. * y * ozi;
5020  curl_shape(7,2) = 0.;
5021 }
5022 
5023 const double Nedelec1PyrFiniteElement::tk[8][3] =
5024 {{1,0,0}, {0,1,0}, {1,0,0}, {0,1,0}, {0,0,1}, {-1,0,1}, {-1,-1,1}, {0,-1,1}};
5025 
5026 void Nedelec1PyrFiniteElement::GetLocalInterpolation (
5027  ElementTransformation &Trans, DenseMatrix &I) const
5028 {
5029  int k, j;
5030 #ifdef MFEM_THREAD_SAFE
5031  DenseMatrix vshape(dof, dim);
5032 #endif
5033 
5034 #ifdef MFEM_DEBUG
5035  for (k = 0; k < dof; k++)
5036  {
5037  CalcVShape (Nodes.IntPoint(k), vshape);
5038  for (j = 0; j < dof; j++)
5039  {
5040  double d = ( vshape(j,0)*tk[k][0] + vshape(j,1)*tk[k][1] +
5041  vshape(j,2)*tk[k][2] );
5042  if (j == k) { d -= 1.0; }
5043  if (fabs(d) > 1.0e-12)
5044  {
5045  mfem::err << "Nedelec1PyrFiniteElement::GetLocalInterpolation (...)\n"
5046  " k = " << k << ", j = " << j << ", d = " << d << endl;
5047  mfem_error();
5048  }
5049  }
5050  }
5051 #endif
5052 
5053  IntegrationPoint ip;
5054  ip.x = ip.y = ip.z = 0.0;
5055  Trans.SetIntPoint (&ip);
5056  // Trans must be linear
5057  const DenseMatrix &J = Trans.Jacobian();
5058  double vk[3];
5059  Vector xk (vk, 3);
5060 
5061  for (k = 0; k < dof; k++)
5062  {
5063  Trans.Transform (Nodes.IntPoint (k), xk);
5064  ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
5065  CalcVShape (ip, vshape);
5066  // vk = J tk
5067  vk[0] = J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2];
5068  vk[1] = J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2];
5069  vk[2] = J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2];
5070  for (j = 0; j < dof; j++)
5071  if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
5072  vshape(j,2)*vk[2])) < 1.0e-12)
5073  {
5074  I(k,j) = 0.0;
5075  }
5076  }
5077 }
5078 
5079 void Nedelec1PyrFiniteElement::Project (
5080  VectorCoefficient &vc, ElementTransformation &Trans,
5081  Vector &dofs) const
5082 {
5083  double vk[3];
5084  Vector xk (vk, 3);
5085 
5086  for (int k = 0; k < dof; k++)
5087  {
5088  Trans.SetIntPoint (&Nodes.IntPoint (k));
5089  const DenseMatrix &J = Trans.Jacobian();
5090 
5091  vc.Eval (xk, Trans, Nodes.IntPoint (k));
5092  // xk^t J tk
5093  dofs(k) =
5094  vk[0] * ( J(0,0)*tk[k][0]+J(0,1)*tk[k][1]+J(0,2)*tk[k][2] ) +
5095  vk[1] * ( J(1,0)*tk[k][0]+J(1,1)*tk[k][1]+J(1,2)*tk[k][2] ) +
5096  vk[2] * ( J(2,0)*tk[k][0]+J(2,1)*tk[k][1]+J(2,2)*tk[k][2] );
5097  }
5098 }
5099 
5100 void Nedelec1PyrFiniteElement::ProjectGrad(const FiniteElement &fe,
5101  ElementTransformation &Trans,
5102  DenseMatrix &grad) const
5103 {
5104  DenseMatrix dshape(fe.GetDof(), 3);
5105  Vector grad_k(fe.GetDof());
5106 
5107  grad.SetSize(dof, fe.GetDof());
5108  for (int k = 0; k < dof; k++)
5109  {
5110  fe.CalcDShape(Nodes.IntPoint(k), dshape);
5111  dshape.Mult(tk[k], grad_k);
5112  for (int j = 0; j < grad_k.Size(); j++)
5113  {
5114  grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
5115  }
5116  }
5117 }
5118 
5119 
5120 RT0HexFiniteElement::RT0HexFiniteElement()
5121  : VectorFiniteElement(3, Geometry::CUBE, 6, 1, H_DIV, FunctionSpace::Qk)
5122 {
5123  // not real nodes ...
5124  // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
5125  Nodes.IntPoint(0).x = 0.5;
5126  Nodes.IntPoint(0).y = 0.5;
5127  Nodes.IntPoint(0).z = 0.0;
5128 
5129  Nodes.IntPoint(1).x = 0.5;
5130  Nodes.IntPoint(1).y = 0.0;
5131  Nodes.IntPoint(1).z = 0.5;
5132 
5133  Nodes.IntPoint(2).x = 1.0;
5134  Nodes.IntPoint(2).y = 0.5;
5135  Nodes.IntPoint(2).z = 0.5;
5136 
5137  Nodes.IntPoint(3).x = 0.5;
5138  Nodes.IntPoint(3).y = 1.0;
5139  Nodes.IntPoint(3).z = 0.5;
5140 
5141  Nodes.IntPoint(4).x = 0.0;
5142  Nodes.IntPoint(4).y = 0.5;
5143  Nodes.IntPoint(4).z = 0.5;
5144 
5145  Nodes.IntPoint(5).x = 0.5;
5146  Nodes.IntPoint(5).y = 0.5;
5147  Nodes.IntPoint(5).z = 1.0;
5148 }
5149 
5150 void RT0HexFiniteElement::CalcVShape(const IntegrationPoint &ip,
5151  DenseMatrix &shape) const
5152 {
5153  double x = ip.x, y = ip.y, z = ip.z;
5154  // z = 0
5155  shape(0,0) = 0.;
5156  shape(0,1) = 0.;
5157  shape(0,2) = z - 1.;
5158  // y = 0
5159  shape(1,0) = 0.;
5160  shape(1,1) = y - 1.;
5161  shape(1,2) = 0.;
5162  // x = 1
5163  shape(2,0) = x;
5164  shape(2,1) = 0.;
5165  shape(2,2) = 0.;
5166  // y = 1
5167  shape(3,0) = 0.;
5168  shape(3,1) = y;
5169  shape(3,2) = 0.;
5170  // x = 0
5171  shape(4,0) = x - 1.;
5172  shape(4,1) = 0.;
5173  shape(4,2) = 0.;
5174  // z = 1
5175  shape(5,0) = 0.;
5176  shape(5,1) = 0.;
5177  shape(5,2) = z;
5178 }
5179 
5180 void RT0HexFiniteElement::CalcDivShape(const IntegrationPoint &ip,
5181  Vector &divshape) const
5182 {
5183  divshape(0) = 1.;
5184  divshape(1) = 1.;
5185  divshape(2) = 1.;
5186  divshape(3) = 1.;
5187  divshape(4) = 1.;
5188  divshape(5) = 1.;
5189 }
5190 
5191 const double RT0HexFiniteElement::nk[6][3] =
5192 {{0,0,-1}, {0,-1,0}, {1,0,0}, {0,1,0}, {-1,0,0}, {0,0,1}};
5193 
5194 void RT0HexFiniteElement::GetLocalInterpolation (
5195  ElementTransformation &Trans, DenseMatrix &I) const
5196 {
5197  int k, j;
5198 #ifdef MFEM_THREAD_SAFE
5199  DenseMatrix vshape(dof, dim);
5200 #endif
5201 
5202 #ifdef MFEM_DEBUG
5203  for (k = 0; k < 6; k++)
5204  {
5205  CalcVShape (Nodes.IntPoint(k), vshape);
5206  for (j = 0; j < 6; j++)
5207  {
5208  double d = ( vshape(j,0)*nk[k][0] + vshape(j,1)*nk[k][1] +
5209  vshape(j,2)*nk[k][2] );
5210  if (j == k) { d -= 1.0; }
5211  if (fabs(d) > 1.0e-12)
5212  {
5213  mfem::err << "RT0HexFiniteElement::GetLocalInterpolation (...)\n"
5214  " k = " << k << ", j = " << j << ", d = " << d << endl;
5215  mfem_error();
5216  }
5217  }
5218  }
5219 #endif
5220 
5221  IntegrationPoint ip;
5222  ip.x = ip.y = ip.z = 0.0;
5223  Trans.SetIntPoint (&ip);
5224  // Trans must be linear
5225  // set Jinv = |J| J^{-t} = adj(J)^t
5226  const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
5227 
5228  double vk[3];
5229  Vector xk (vk, 3);
5230 
5231  for (k = 0; k < 6; k++)
5232  {
5233  Trans.Transform (Nodes.IntPoint (k), xk);
5234  ip.x = vk[0]; ip.y = vk[1]; ip.z = vk[2];
5235  CalcVShape (ip, vshape);
5236  // vk = |J| J^{-t} nk
5237  vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2];
5238  vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2];
5239  vk[2] = Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2];
5240  for (j = 0; j < 6; j++)
5241  if (fabs (I(k,j) = (vshape(j,0)*vk[0]+vshape(j,1)*vk[1]+
5242  vshape(j,2)*vk[2])) < 1.0e-12)
5243  {
5244  I(k,j) = 0.0;
5245  }
5246  }
5247 }
5248 
5249 void RT0HexFiniteElement::Project (
5250  VectorCoefficient &vc, ElementTransformation &Trans,
5251  Vector &dofs) const
5252 {
5253  double vk[3];
5254  Vector xk (vk, 3);
5255 
5256  for (int k = 0; k < 6; k++)
5257  {
5258  Trans.SetIntPoint (&Nodes.IntPoint (k));
5259  // set Jinv = |J| J^{-t} = adj(J)^t
5260  const DenseMatrix &Jinv = Trans.TransposeAdjugateJacobian();
5261 
5262  vc.Eval (xk, Trans, Nodes.IntPoint (k));
5263  // xk^t |J| J^{-t} nk
5264  dofs(k) =
5265  vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1]+Jinv(0,2)*nk[k][2] ) +
5266  vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1]+Jinv(1,2)*nk[k][2] ) +
5267  vk[2] * ( Jinv(2,0)*nk[k][0]+Jinv(2,1)*nk[k][1]+Jinv(2,2)*nk[k][2] );
5268  }
5269 }
5270 
5271 RT1HexFiniteElement::RT1HexFiniteElement()
5272  : VectorFiniteElement(3, Geometry::CUBE, 36, 2, H_DIV,
5273  FunctionSpace::Qk)
5274 {
5275  // z = 0
5276  Nodes.IntPoint(2).x = 1./3.;
5277  Nodes.IntPoint(2).y = 1./3.;
5278  Nodes.IntPoint(2).z = 0.0;
5279  Nodes.IntPoint(3).x = 2./3.;
5280  Nodes.IntPoint(3).y = 1./3.;
5281  Nodes.IntPoint(3).z = 0.0;
5282  Nodes.IntPoint(0).x = 1./3.;
5283  Nodes.IntPoint(0).y = 2./3.;
5284  Nodes.IntPoint(0).z = 0.0;
5285  Nodes.IntPoint(1).x = 2./3.;
5286  Nodes.IntPoint(1).y = 2./3.;
5287  Nodes.IntPoint(1).z = 0.0;
5288  // y = 0
5289  Nodes.IntPoint(4).x = 1./3.;
5290  Nodes.IntPoint(4).y = 0.0;
5291  Nodes.IntPoint(4).z = 1./3.;
5292  Nodes.IntPoint(5).x = 2./3.;
5293  Nodes.IntPoint(5).y = 0.0;
5294  Nodes.IntPoint(5).z = 1./3.;
5295  Nodes.IntPoint(6).x = 1./3.;
5296  Nodes.IntPoint(6).y = 0.0;
5297  Nodes.IntPoint(6).z = 2./3.;
5298  Nodes.IntPoint(7).x = 2./3.;
5299  Nodes.IntPoint(7).y = 0.0;
5300  Nodes.IntPoint(7).z = 2./3.;
5301  // x = 1
5302  Nodes.IntPoint(8).x = 1.0;
5303  Nodes.IntPoint(8).y = 1./3.;
5304  Nodes.IntPoint(8).z = 1./3.;
5305  Nodes.IntPoint(9).x = 1.0;
5306  Nodes.IntPoint(9).y = 2./3.;
5307  Nodes.IntPoint(9).z = 1./3.;
5308  Nodes.IntPoint(10).x = 1.0;
5309  Nodes.IntPoint(10).y = 1./3.;
5310  Nodes.IntPoint(10).z = 2./3.;
5311  Nodes.IntPoint(11).x = 1.0;
5312  Nodes.IntPoint(11).y = 2./3.;
5313  Nodes.IntPoint(11).z = 2./3.;
5314  // y = 1
5315  Nodes.IntPoint(13).x = 1./3.;
5316  Nodes.IntPoint(13).y = 1.0;
5317  Nodes.IntPoint(13).z = 1./3.;
5318  Nodes.IntPoint(12).x = 2./3.;
5319  Nodes.IntPoint(12).y = 1.0;
5320  Nodes.IntPoint(12).z = 1./3.;
5321  Nodes.IntPoint(15).x = 1./3.;
5322  Nodes.IntPoint(15).y = 1.0;
5323  Nodes.IntPoint(15).z = 2./3.;
5324  Nodes.IntPoint(14).x = 2./3.;
5325  Nodes.IntPoint(14).y = 1.0;
5326  Nodes.IntPoint(14).z = 2./3.;
5327  // x = 0
5328  Nodes.IntPoint(17).x = 0.0;
5329  Nodes.IntPoint(17).y = 1./3.;
5330  Nodes.IntPoint(17).z = 1./3.;
5331  Nodes.IntPoint(16).x = 0.0;
5332  Nodes.IntPoint(16).y = 2./3.;
5333  Nodes.IntPoint(16).z = 1./3.;
5334  Nodes.IntPoint(19).x = 0.0;
5335  Nodes.IntPoint(19).y = 1./3.;
5336  Nodes.IntPoint(19).z = 2./3.;
5337  Nodes.IntPoint(18).x = 0.0;
5338  Nodes.IntPoint(18).y = 2./3.;
5339  Nodes.IntPoint(18).z = 2./3.;
5340  // z = 1
5341  Nodes.IntPoint(20).x = 1./3.;
5342  Nodes.IntPoint(20).y = 1./3.;
5343  Nodes.IntPoint(20).z = 1.0;
5344  Nodes.IntPoint(21).x = 2./3.;
5345  Nodes.IntPoint(21).y = 1./3.;
5346  Nodes.IntPoint(21).z = 1.0;
5347  Nodes.IntPoint(22).x = 1./3.;
5348  Nodes.IntPoint(22).y = 2./3.;
5349  Nodes.IntPoint(22).z = 1.0;
5350  Nodes.IntPoint(23).x = 2./3.;
5351  Nodes.IntPoint(23).y = 2./3.;
5352  Nodes.IntPoint(23).z = 1.0;
5353  // x = 0.5 (interior)
5354  Nodes.IntPoint(24).x = 0.5;
5355  Nodes.IntPoint(24).y = 1./3.;
5356  Nodes.IntPoint(24).z = 1./3.;
5357  Nodes.IntPoint(25).x = 0.5;
5358  Nodes.IntPoint(25).y = 1./3.;
5359  Nodes.IntPoint(25).z = 2./3.;
5360  Nodes.IntPoint(26).x = 0.5;
5361  Nodes.IntPoint(26).y = 2./3.;
5362  Nodes.IntPoint(26).z = 1./3.;
5363  Nodes.IntPoint(27).x = 0.5;
5364  Nodes.IntPoint(27).y = 2./3.;
5365  Nodes.IntPoint(27).z = 2./3.;
5366  // y = 0.5 (interior)
5367  Nodes.IntPoint(28).x = 1./3.;
5368  Nodes.IntPoint(28).y = 0.5;
5369  Nodes.IntPoint(28).z = 1./3.;
5370  Nodes.IntPoint(29).x = 1./3.;
5371  Nodes.IntPoint(29).y = 0.5;
5372  Nodes.IntPoint(29).z = 2./3.;
5373  Nodes.IntPoint(30).x = 2./3.;
5374  Nodes.IntPoint(30).y = 0.5;
5375  Nodes.IntPoint(30).z = 1./3.;
5376  Nodes.IntPoint(31).x = 2./3.;
5377  Nodes.IntPoint(31).y = 0.5;
5378  Nodes.IntPoint(31).z = 2./3.;
5379  // z = 0.5 (interior)
5380  Nodes.IntPoint(32).x = 1./3.;
5381  Nodes.IntPoint(32).y = 1./3.;
5382  Nodes.IntPoint(32).z = 0.5;
5383  Nodes.IntPoint(33).x = 1./3.;
5384  Nodes.IntPoint(33).y = 2./3.;
5385  Nodes.IntPoint(33).z = 0.5;
5386  Nodes.IntPoint(34).x = 2./3.;
5387  Nodes.IntPoint(34).y = 1./3.;
5388  Nodes.IntPoint(34).z = 0.5;
5389  Nodes.IntPoint(35).x = 2./3.;
5390  Nodes.IntPoint(35).y = 2./3.;
5391  Nodes.IntPoint(35).z = 0.5;
5392 }
5393 
5394 void RT1HexFiniteElement::CalcVShape(const IntegrationPoint &ip,
5395  DenseMatrix &shape) const
5396 {
5397  double x = ip.x, y = ip.y, z = ip.z;
5398  // z = 0
5399  shape(2,0) = 0.;
5400  shape(2,1) = 0.;
5401  shape(2,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5402  shape(3,0) = 0.;
5403  shape(3,1) = 0.;
5404  shape(3,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5405  shape(0,0) = 0.;
5406  shape(0,1) = 0.;
5407  shape(0,2) = -(1. - 3.*z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5408  shape(1,0) = 0.;
5409  shape(1,1) = 0.;
5410  shape(1,2) = -(1. - 3.*z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5411  // y = 0
5412  shape(4,0) = 0.;
5413  shape(4,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5414  shape(4,2) = 0.;
5415  shape(5,0) = 0.;
5416  shape(5,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5417  shape(5,2) = 0.;
5418  shape(6,0) = 0.;
5419  shape(6,1) = -(1. - 3.*y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5420  shape(6,2) = 0.;
5421  shape(7,0) = 0.;
5422  shape(7,1) = -(1. - 3.*y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5423  shape(7,2) = 0.;
5424  // x = 1
5425  shape(8,0) = (-x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5426  shape(8,1) = 0.;
5427  shape(8,2) = 0.;
5428  shape(9,0) = (-x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5429  shape(9,1) = 0.;
5430  shape(9,2) = 0.;
5431  shape(10,0) = (-x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5432  shape(10,1) = 0.;
5433  shape(10,2) = 0.;
5434  shape(11,0) = (-x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5435  shape(11,1) = 0.;
5436  shape(11,2) = 0.;
5437  // y = 1
5438  shape(13,0) = 0.;
5439  shape(13,1) = (-y + 2.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5440  shape(13,2) = 0.;
5441  shape(12,0) = 0.;
5442  shape(12,1) = (-y + 2.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5443  shape(12,2) = 0.;
5444  shape(15,0) = 0.;
5445  shape(15,1) = (-y + 2.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5446  shape(15,2) = 0.;
5447  shape(14,0) = 0.;
5448  shape(14,1) = (-y + 2.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5449  shape(14,2) = 0.;
5450  // x = 0
5451  shape(17,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5452  shape(17,1) = 0.;
5453  shape(17,2) = 0.;
5454  shape(16,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5455  shape(16,1) = 0.;
5456  shape(16,2) = 0.;
5457  shape(19,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5458  shape(19,1) = 0.;
5459  shape(19,2) = 0.;
5460  shape(18,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5461  shape(18,1) = 0.;
5462  shape(18,2) = 0.;
5463  // z = 1
5464  shape(20,0) = 0.;
5465  shape(20,1) = 0.;
5466  shape(20,2) = (-z + 2.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5467  shape(21,0) = 0.;
5468  shape(21,1) = 0.;
5469  shape(21,2) = (-z + 2.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5470  shape(22,0) = 0.;
5471  shape(22,1) = 0.;
5472  shape(22,2) = (-z + 2.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5473  shape(23,0) = 0.;
5474  shape(23,1) = 0.;
5475  shape(23,2) = (-z + 2.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5476  // x = 0.5 (interior)
5477  shape(24,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*( 2. - 3.*z);
5478  shape(24,1) = 0.;
5479  shape(24,2) = 0.;
5480  shape(25,0) = (4.*x - 4.*x*x)*( 2. - 3.*y)*(-1. + 3.*z);
5481  shape(25,1) = 0.;
5482  shape(25,2) = 0.;
5483  shape(26,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*( 2. - 3.*z);
5484  shape(26,1) = 0.;
5485  shape(26,2) = 0.;
5486  shape(27,0) = (4.*x - 4.*x*x)*(-1. + 3.*y)*(-1. + 3.*z);
5487  shape(27,1) = 0.;
5488  shape(27,2) = 0.;
5489  // y = 0.5 (interior)
5490  shape(28,0) = 0.;
5491  shape(28,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*( 2. - 3.*z);
5492  shape(28,2) = 0.;
5493  shape(29,0) = 0.;
5494  shape(29,1) = (4.*y - 4.*y*y)*( 2. - 3.*x)*(-1. + 3.*z);
5495  shape(29,2) = 0.;
5496  shape(30,0) = 0.;
5497  shape(30,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*( 2. - 3.*z);
5498  shape(30,2) = 0.;
5499  shape(31,0) = 0.;
5500  shape(31,1) = (4.*y - 4.*y*y)*(-1. + 3.*x)*(-1. + 3.*z);
5501  shape(31,2) = 0.;
5502  // z = 0.5 (interior)
5503  shape(32,0) = 0.;
5504  shape(32,1) = 0.;
5505  shape(32,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*( 2. - 3.*y);
5506  shape(33,0) = 0.;
5507  shape(33,1) = 0.;
5508  shape(33,2) = (4.*z - 4.*z*z)*( 2. - 3.*x)*(-1. + 3.*y);
5509  shape(34,0) = 0.;
5510  shape(34,1) = 0.;
5511  shape(34,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*( 2. - 3.*y);
5512  shape(35,0) = 0.;
5513  shape(35,1) = 0.;
5514  shape(35,2) = (4.*z - 4.*z*z)*(-1. + 3.*x)*(-1. + 3.*y);
5515 }
5516 
5517 void RT1HexFiniteElement::CalcDivShape(const IntegrationPoint &ip,
5518  Vector &divshape) const
5519 {
5520  double x = ip.x, y = ip.y, z = ip.z;
5521  // z = 0
5522  divshape(2) = -(-3. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5523  divshape(3) = -(-3. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5524  divshape(0) = -(-3. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5525  divshape(1) = -(-3. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5526  // y = 0
5527  divshape(4) = -(-3. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5528  divshape(5) = -(-3. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5529  divshape(6) = -(-3. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5530  divshape(7) = -(-3. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5531  // x = 1
5532  divshape(8) = (-1. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5533  divshape(9) = (-1. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5534  divshape(10) = (-1. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5535  divshape(11) = (-1. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5536  // y = 1
5537  divshape(13) = (-1. + 4.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5538  divshape(12) = (-1. + 4.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5539  divshape(15) = (-1. + 4.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5540  divshape(14) = (-1. + 4.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5541  // x = 0
5542  divshape(17) = -(-3. + 4.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5543  divshape(16) = -(-3. + 4.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5544  divshape(19) = -(-3. + 4.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5545  divshape(18) = -(-3. + 4.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5546  // z = 1
5547  divshape(20) = (-1. + 4.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5548  divshape(21) = (-1. + 4.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5549  divshape(22) = (-1. + 4.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5550  divshape(23) = (-1. + 4.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5551  // x = 0.5 (interior)
5552  divshape(24) = ( 4. - 8.*x)*( 2. - 3.*y)*( 2. - 3.*z);
5553  divshape(25) = ( 4. - 8.*x)*( 2. - 3.*y)*(-1. + 3.*z);
5554  divshape(26) = ( 4. - 8.*x)*(-1. + 3.*y)*( 2. - 3.*z);
5555  divshape(27) = ( 4. - 8.*x)*(-1. + 3.*y)*(-1. + 3.*z);
5556  // y = 0.5 (interior)
5557  divshape(28) = ( 4. - 8.*y)*( 2. - 3.*x)*( 2. - 3.*z);
5558  divshape(29) = ( 4. - 8.*y)*( 2. - 3.*x)*(-1. + 3.*z);
5559  divshape(30) = ( 4. - 8.*y)*(-1. + 3.*x)*( 2. - 3.*z);
5560  divshape(31) = ( 4. - 8.*y)*(-1. + 3.*x)*(-1. + 3.*z);
5561  // z = 0.5 (interior)
5562  divshape(32) = ( 4. - 8.*z)*( 2. - 3.*x)*( 2. - 3.*y);
5563  divshape(33) = ( 4. - 8.*z)*( 2. - 3.*x)*(-1. + 3.*y);
5564  divshape(34) = ( 4. - 8.*z)*(-1. + 3.*x)*( 2. - 3.*y);
5565  divshape(35) = ( 4. - 8.*z)*(-1. + 3.*x)*(-1. + 3.*y);
5566 }
5567 
5568 const double RT1HexFiniteElement::nk[36][3] =
5569 {
5570  {0, 0,-1}, {0, 0,-1}, {0, 0,-1}, {0, 0,-1},
5571  {0,-1, 0}, {0,-1, 0}, {0,-1, 0}, {0,-1, 0},
5572  {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5573  {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5574  {-1,0, 0}, {-1,0, 0}, {-1,0, 0}, {-1,0, 0},
5575  {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1},
5576  {1, 0, 0}, {1, 0, 0}, {1, 0, 0}, {1, 0, 0},
5577  {0, 1, 0}, {0, 1, 0}, {0, 1, 0}, {0, 1, 0},
5578  {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}
5579 };
5580 
5581 void RT1HexFiniteElement::GetLocalInterpolation (
5582  ElementTransformation &Trans, DenseMatrix &I) const
5583 {
5584  int k, j;
5585 #ifdef MFEM_THREAD_SAFE
5586  DenseMatrix vshape(dof, dim);
5587 #endif
5588 
5589 #ifdef MFEM_DEBUG
5590  for (k = 0; k < 36; k++)
5591  {
5592  CalcVShape (Nodes.IntPoint(k), vshape);
5593  for (j = 0; j < 36; j++)
5594  {
5595  double d = ( vshape(j,0)*nk[k][0] + vshape(j,1)*nk[k][1] +
5596  vshape(j,2)*nk[k][2] );
5597  if (j == k) { d -= 1.0; }
5598  if (fabs(d) > 1.0e-12)
5599  {
5600  mfem::err << "RT0HexFiniteElement::GetLocal