MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fe_nd.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 // Nedelec Finite Element classes
13 
14 #include "fe_nd.hpp"
15 #include "../coefficient.hpp"
16 
17 namespace mfem
18 {
19 
20 using namespace std;
21 
22 const double ND_HexahedronElement::tk[18] =
23 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,0.,0., 0.,-1.,0., 0.,0.,-1. };
24 
26  const int cb_type, const int ob_type)
27  : VectorTensorFiniteElement(3, 3*p*(p + 1)*(p + 1), p, cb_type, ob_type,
28  H_CURL, DofMapType::L2_DOF_MAP),
29  dof2tk(dof), cp(poly1d.ClosedPoints(p, cb_type))
30 {
31  if (obasis1d.IsIntegratedType()) { is_nodal = false; }
32 
34 
35  const double *op = poly1d.OpenPoints(p - 1, ob_type);
36  const int dof3 = dof/3;
37 
38 #ifndef MFEM_THREAD_SAFE
39  shape_cx.SetSize(p + 1);
40  shape_ox.SetSize(p);
41  shape_cy.SetSize(p + 1);
42  shape_oy.SetSize(p);
43  shape_cz.SetSize(p + 1);
44  shape_oz.SetSize(p);
45  dshape_cx.SetSize(p + 1);
46  dshape_cy.SetSize(p + 1);
47  dshape_cz.SetSize(p + 1);
48 #endif
49 
50  // edges
51  int o = 0;
52  for (int i = 0; i < p; i++) // (0,1)
53  {
54  dof_map[0*dof3 + i + (0 + 0*(p + 1))*p] = o++;
55  }
56  for (int i = 0; i < p; i++) // (1,2)
57  {
58  dof_map[1*dof3 + p + (i + 0*p)*(p + 1)] = o++;
59  }
60  for (int i = 0; i < p; i++) // (3,2)
61  {
62  dof_map[0*dof3 + i + (p + 0*(p + 1))*p] = o++;
63  }
64  for (int i = 0; i < p; i++) // (0,3)
65  {
66  dof_map[1*dof3 + 0 + (i + 0*p)*(p + 1)] = o++;
67  }
68  for (int i = 0; i < p; i++) // (4,5)
69  {
70  dof_map[0*dof3 + i + (0 + p*(p + 1))*p] = o++;
71  }
72  for (int i = 0; i < p; i++) // (5,6)
73  {
74  dof_map[1*dof3 + p + (i + p*p)*(p + 1)] = o++;
75  }
76  for (int i = 0; i < p; i++) // (7,6)
77  {
78  dof_map[0*dof3 + i + (p + p*(p + 1))*p] = o++;
79  }
80  for (int i = 0; i < p; i++) // (4,7)
81  {
82  dof_map[1*dof3 + 0 + (i + p*p)*(p + 1)] = o++;
83  }
84  for (int i = 0; i < p; i++) // (0,4)
85  {
86  dof_map[2*dof3 + 0 + (0 + i*(p + 1))*(p + 1)] = o++;
87  }
88  for (int i = 0; i < p; i++) // (1,5)
89  {
90  dof_map[2*dof3 + p + (0 + i*(p + 1))*(p + 1)] = o++;
91  }
92  for (int i = 0; i < p; i++) // (2,6)
93  {
94  dof_map[2*dof3 + p + (p + i*(p + 1))*(p + 1)] = o++;
95  }
96  for (int i = 0; i < p; i++) // (3,7)
97  {
98  dof_map[2*dof3 + 0 + (p + i*(p + 1))*(p + 1)] = o++;
99  }
100 
101  // faces
102  // (3,2,1,0) -- bottom
103  for (int j = 1; j < p; j++) // x - components
104  for (int i = 0; i < p; i++)
105  {
106  dof_map[0*dof3 + i + ((p - j) + 0*(p + 1))*p] = o++;
107  }
108  for (int j = 0; j < p; j++) // y - components
109  for (int i = 1; i < p; i++)
110  {
111  dof_map[1*dof3 + i + ((p - 1 - j) + 0*p)*(p + 1)] = -1 - (o++);
112  }
113  // (0,1,5,4) -- front
114  for (int k = 1; k < p; k++) // x - components
115  for (int i = 0; i < p; i++)
116  {
117  dof_map[0*dof3 + i + (0 + k*(p + 1))*p] = o++;
118  }
119  for (int k = 0; k < p; k++) // z - components
120  for (int i = 1; i < p; i++ )
121  {
122  dof_map[2*dof3 + i + (0 + k*(p + 1))*(p + 1)] = o++;
123  }
124  // (1,2,6,5) -- right
125  for (int k = 1; k < p; k++) // y - components
126  for (int j = 0; j < p; j++)
127  {
128  dof_map[1*dof3 + p + (j + k*p)*(p + 1)] = o++;
129  }
130  for (int k = 0; k < p; k++) // z - components
131  for (int j = 1; j < p; j++)
132  {
133  dof_map[2*dof3 + p + (j + k*(p + 1))*(p + 1)] = o++;
134  }
135  // (2,3,7,6) -- back
136  for (int k = 1; k < p; k++) // x - components
137  for (int i = 0; i < p; i++)
138  {
139  dof_map[0*dof3 + (p - 1 - i) + (p + k*(p + 1))*p] = -1 - (o++);
140  }
141  for (int k = 0; k < p; k++) // z - components
142  for (int i = 1; i < p; i++)
143  {
144  dof_map[2*dof3 + (p - i) + (p + k*(p + 1))*(p + 1)] = o++;
145  }
146  // (3,0,4,7) -- left
147  for (int k = 1; k < p; k++) // y - components
148  for (int j = 0; j < p; j++)
149  {
150  dof_map[1*dof3 + 0 + ((p - 1 - j) + k*p)*(p + 1)] = -1 - (o++);
151  }
152  for (int k = 0; k < p; k++) // z - components
153  for (int j = 1; j < p; j++)
154  {
155  dof_map[2*dof3 + 0 + ((p - j) + k*(p + 1))*(p + 1)] = o++;
156  }
157  // (4,5,6,7) -- top
158  for (int j = 1; j < p; j++) // x - components
159  for (int i = 0; i < p; i++)
160  {
161  dof_map[0*dof3 + i + (j + p*(p + 1))*p] = o++;
162  }
163  for (int j = 0; j < p; j++) // y - components
164  for (int i = 1; i < p; i++)
165  {
166  dof_map[1*dof3 + i + (j + p*p)*(p + 1)] = o++;
167  }
168 
169  // interior
170  // x-components
171  for (int k = 1; k < p; k++)
172  for (int j = 1; j < p; j++)
173  for (int i = 0; i < p; i++)
174  {
175  dof_map[0*dof3 + i + (j + k*(p + 1))*p] = o++;
176  }
177  // y-components
178  for (int k = 1; k < p; k++)
179  for (int j = 0; j < p; j++)
180  for (int i = 1; i < p; i++)
181  {
182  dof_map[1*dof3 + i + (j + k*p)*(p + 1)] = o++;
183  }
184  // z-components
185  for (int k = 0; k < p; k++)
186  for (int j = 1; j < p; j++)
187  for (int i = 1; i < p; i++)
188  {
189  dof_map[2*dof3 + i + (j + k*(p + 1))*(p + 1)] = o++;
190  }
191 
192  // set dof2tk and Nodes
193  o = 0;
194  // x-components
195  for (int k = 0; k <= p; k++)
196  for (int j = 0; j <= p; j++)
197  for (int i = 0; i < p; i++)
198  {
199  int idx;
200  if ((idx = dof_map[o++]) < 0)
201  {
202  dof2tk[idx = -1 - idx] = 3;
203  }
204  else
205  {
206  dof2tk[idx] = 0;
207  }
208  Nodes.IntPoint(idx).Set3(op[i], cp[j], cp[k]);
209  }
210  // y-components
211  for (int k = 0; k <= p; k++)
212  for (int j = 0; j < p; j++)
213  for (int i = 0; i <= p; i++)
214  {
215  int idx;
216  if ((idx = dof_map[o++]) < 0)
217  {
218  dof2tk[idx = -1 - idx] = 4;
219  }
220  else
221  {
222  dof2tk[idx] = 1;
223  }
224  Nodes.IntPoint(idx).Set3(cp[i], op[j], cp[k]);
225  }
226  // z-components
227  for (int k = 0; k < p; k++)
228  for (int j = 0; j <= p; j++)
229  for (int i = 0; i <= p; i++)
230  {
231  int idx;
232  if ((idx = dof_map[o++]) < 0)
233  {
234  dof2tk[idx = -1 - idx] = 5;
235  }
236  else
237  {
238  dof2tk[idx] = 2;
239  }
240  Nodes.IntPoint(idx).Set3(cp[i], cp[j], op[k]);
241  }
242 }
243 
246  Vector &dofs) const
247 {
248  MFEM_ASSERT(obasis1d.IsIntegratedType(), "Not integrated type");
249  double vk[Geometry::MaxDim];
250  Vector xk(vk, vc.GetVDim());
251 
253  const int nqpt = ir.GetNPoints();
254 
255  IntegrationPoint ip3d;
256 
257  int o = 0;
258  for (int c = 0; c < 3; ++c) // loop over x, y, z components
259  {
260  const int im = c == 0 ? order - 1 : order;
261  const int jm = c == 1 ? order - 1 : order;
262  const int km = c == 2 ? order - 1 : order;
263 
264  for (int k = 0; k <= km; k++)
265  for (int j = 0; j <= jm; j++)
266  for (int i = 0; i <= im; i++)
267  {
268  int idx;
269  if ((idx = dof_map[o++]) < 0)
270  {
271  idx = -1 - idx;
272  }
273 
274  const int id1 = c == 0 ? i : (c == 1 ? j : k);
275  const double h = cp[id1+1] - cp[id1];
276 
277  double val = 0.0;
278 
279  for (int q = 0; q < nqpt; q++)
280  {
281  const IntegrationPoint &ip1d = ir.IntPoint(q);
282 
283  if (c == 0)
284  {
285  ip3d.Set3(cp[i] + (h*ip1d.x), cp[j], cp[k]);
286  }
287  else if (c == 1)
288  {
289  ip3d.Set3(cp[i], cp[j] + (h*ip1d.x), cp[k]);
290  }
291  else
292  {
293  ip3d.Set3(cp[i], cp[j], cp[k] + (h*ip1d.x));
294  }
295 
296  Trans.SetIntPoint(&ip3d);
297  vc.Eval(xk, Trans, ip3d);
298 
299  // xk^t J tk
300  const double ipval = Trans.Jacobian().InnerProduct(tk + dof2tk[idx]*dim, vk);
301  val += ip1d.weight * ipval;
302  }
303 
304  dofs(idx) = val*h;
305  }
306  }
307 }
308 
310  DenseMatrix &shape) const
311 {
312  const int p = order;
313 
314 #ifdef MFEM_THREAD_SAFE
315  Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
316  Vector shape_cz(p + 1), shape_oz(p);
317  Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
318 #endif
319 
321  {
322  cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
323  cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
324  cbasis1d.Eval(ip.z, shape_cz, dshape_cz);
325  obasis1d.ScaleIntegrated(false);
326  obasis1d.EvalIntegrated(dshape_cx, shape_ox);
327  obasis1d.EvalIntegrated(dshape_cy, shape_oy);
328  obasis1d.EvalIntegrated(dshape_cz, shape_oz);
329  }
330  else
331  {
332  cbasis1d.Eval(ip.x, shape_cx);
333  cbasis1d.Eval(ip.y, shape_cy);
334  cbasis1d.Eval(ip.z, shape_cz);
335  obasis1d.Eval(ip.x, shape_ox);
336  obasis1d.Eval(ip.y, shape_oy);
337  obasis1d.Eval(ip.z, shape_oz);
338  }
339 
340  int o = 0;
341  // x-components
342  for (int k = 0; k <= p; k++)
343  for (int j = 0; j <= p; j++)
344  for (int i = 0; i < p; i++)
345  {
346  int idx, s;
347  if ((idx = dof_map[o++]) < 0)
348  {
349  idx = -1 - idx, s = -1;
350  }
351  else
352  {
353  s = +1;
354  }
355  shape(idx,0) = s*shape_ox(i)*shape_cy(j)*shape_cz(k);
356  shape(idx,1) = 0.;
357  shape(idx,2) = 0.;
358  }
359  // y-components
360  for (int k = 0; k <= p; k++)
361  for (int j = 0; j < p; j++)
362  for (int i = 0; i <= p; i++)
363  {
364  int idx, s;
365  if ((idx = dof_map[o++]) < 0)
366  {
367  idx = -1 - idx, s = -1;
368  }
369  else
370  {
371  s = +1;
372  }
373  shape(idx,0) = 0.;
374  shape(idx,1) = s*shape_cx(i)*shape_oy(j)*shape_cz(k);
375  shape(idx,2) = 0.;
376  }
377  // z-components
378  for (int k = 0; k < p; k++)
379  for (int j = 0; j <= p; j++)
380  for (int i = 0; i <= p; i++)
381  {
382  int idx, s;
383  if ((idx = dof_map[o++]) < 0)
384  {
385  idx = -1 - idx, s = -1;
386  }
387  else
388  {
389  s = +1;
390  }
391  shape(idx,0) = 0.;
392  shape(idx,1) = 0.;
393  shape(idx,2) = s*shape_cx(i)*shape_cy(j)*shape_oz(k);
394  }
395 }
396 
398  DenseMatrix &curl_shape) const
399 {
400  const int p = order;
401 
402 #ifdef MFEM_THREAD_SAFE
403  Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
404  Vector shape_cz(p + 1), shape_oz(p);
405  Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
406 #endif
407 
408  cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
409  cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
410  cbasis1d.Eval(ip.z, shape_cz, dshape_cz);
412  {
413  obasis1d.ScaleIntegrated(false);
414  obasis1d.EvalIntegrated(dshape_cx, shape_ox);
415  obasis1d.EvalIntegrated(dshape_cy, shape_oy);
416  obasis1d.EvalIntegrated(dshape_cz, shape_oz);
417  }
418  else
419  {
420  obasis1d.Eval(ip.x, shape_ox);
421  obasis1d.Eval(ip.y, shape_oy);
422  obasis1d.Eval(ip.z, shape_oz);
423  }
424 
425  int o = 0;
426  // x-components
427  for (int k = 0; k <= p; k++)
428  for (int j = 0; j <= p; j++)
429  for (int i = 0; i < p; i++)
430  {
431  int idx, s;
432  if ((idx = dof_map[o++]) < 0)
433  {
434  idx = -1 - idx, s = -1;
435  }
436  else
437  {
438  s = +1;
439  }
440  curl_shape(idx,0) = 0.;
441  curl_shape(idx,1) = s*shape_ox(i)* shape_cy(j)*dshape_cz(k);
442  curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j)* shape_cz(k);
443  }
444  // y-components
445  for (int k = 0; k <= p; k++)
446  for (int j = 0; j < p; j++)
447  for (int i = 0; i <= p; i++)
448  {
449  int idx, s;
450  if ((idx = dof_map[o++]) < 0)
451  {
452  idx = -1 - idx, s = -1;
453  }
454  else
455  {
456  s = +1;
457  }
458  curl_shape(idx,0) = -s* shape_cx(i)*shape_oy(j)*dshape_cz(k);
459  curl_shape(idx,1) = 0.;
460  curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j)* shape_cz(k);
461  }
462  // z-components
463  for (int k = 0; k < p; k++)
464  for (int j = 0; j <= p; j++)
465  for (int i = 0; i <= p; i++)
466  {
467  int idx, s;
468  if ((idx = dof_map[o++]) < 0)
469  {
470  idx = -1 - idx, s = -1;
471  }
472  else
473  {
474  s = +1;
475  }
476  curl_shape(idx,0) = s* shape_cx(i)*dshape_cy(j)*shape_oz(k);
477  curl_shape(idx,1) = -s*dshape_cx(i)* shape_cy(j)*shape_oz(k);
478  curl_shape(idx,2) = 0.;
479  }
480 }
481 
482 const double ND_QuadrilateralElement::tk[8] =
483 { 1.,0., 0.,1., -1.,0., 0.,-1. };
484 
486  const int cb_type,
487  const int ob_type)
488  : VectorTensorFiniteElement(2, 2*p*(p + 1), p, cb_type, ob_type,
489  H_CURL, DofMapType::L2_DOF_MAP),
490  dof2tk(dof),
491  cp(poly1d.ClosedPoints(p, cb_type))
492 {
493  if (obasis1d.IsIntegratedType()) { is_nodal = false; }
494 
496 
497  const double *op = poly1d.OpenPoints(p - 1, ob_type);
498  const int dof2 = dof/2;
499 
500 #ifndef MFEM_THREAD_SAFE
501  shape_cx.SetSize(p + 1);
502  shape_ox.SetSize(p);
503  shape_cy.SetSize(p + 1);
504  shape_oy.SetSize(p);
505  dshape_cx.SetSize(p + 1);
506  dshape_cy.SetSize(p + 1);
507 #endif
508 
509  // edges
510  int o = 0;
511  for (int i = 0; i < p; i++) // (0,1)
512  {
513  dof_map[0*dof2 + i + 0*p] = o++;
514  }
515  for (int j = 0; j < p; j++) // (1,2)
516  {
517  dof_map[1*dof2 + p + j*(p + 1)] = o++;
518  }
519  for (int i = 0; i < p; i++) // (2,3)
520  {
521  dof_map[0*dof2 + (p - 1 - i) + p*p] = -1 - (o++);
522  }
523  for (int j = 0; j < p; j++) // (3,0)
524  {
525  dof_map[1*dof2 + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
526  }
527 
528  // interior
529  // x-components
530  for (int j = 1; j < p; j++)
531  for (int i = 0; i < p; i++)
532  {
533  dof_map[0*dof2 + i + j*p] = o++;
534  }
535  // y-components
536  for (int j = 0; j < p; j++)
537  for (int i = 1; i < p; i++)
538  {
539  dof_map[1*dof2 + i + j*(p + 1)] = o++;
540  }
541 
542  // set dof2tk and Nodes
543  o = 0;
544  // x-components
545  for (int j = 0; j <= p; j++)
546  for (int i = 0; i < p; i++)
547  {
548  int idx;
549  if ((idx = dof_map[o++]) < 0)
550  {
551  dof2tk[idx = -1 - idx] = 2;
552  }
553  else
554  {
555  dof2tk[idx] = 0;
556  }
557  Nodes.IntPoint(idx).Set2(op[i], cp[j]);
558  }
559  // y-components
560  for (int j = 0; j < p; j++)
561  for (int i = 0; i <= p; i++)
562  {
563  int idx;
564  if ((idx = dof_map[o++]) < 0)
565  {
566  dof2tk[idx = -1 - idx] = 3;
567  }
568  else
569  {
570  dof2tk[idx] = 1;
571  }
572  Nodes.IntPoint(idx).Set2(cp[i], op[j]);
573  }
574 }
575 
578  Vector &dofs) const
579 {
580  MFEM_ASSERT(obasis1d.IsIntegratedType(), "Not integrated type");
581  double vk[Geometry::MaxDim];
582  Vector xk(vk, vc.GetVDim());
583 
585  const int nqpt = ir.GetNPoints();
586 
587  IntegrationPoint ip2d;
588 
589  int o = 0;
590  // x-components
591  for (int j = 0; j <= order; j++)
592  for (int i = 0; i < order; i++)
593  {
594  int idx;
595  if ((idx = dof_map[o++]) < 0)
596  {
597  idx = -1 - idx;
598  }
599 
600  const double h = cp[i+1] - cp[i];
601 
602  double val = 0.0;
603 
604  for (int k = 0; k < nqpt; k++)
605  {
606  const IntegrationPoint &ip1d = ir.IntPoint(k);
607 
608  ip2d.Set2(cp[i] + (h*ip1d.x), cp[j]);
609 
610  Trans.SetIntPoint(&ip2d);
611  vc.Eval(xk, Trans, ip2d);
612 
613  // xk^t J tk
614  const double ipval = Trans.Jacobian().InnerProduct(tk + dof2tk[idx]*dim, vk);
615  val += ip1d.weight * ipval;
616  }
617 
618  dofs(idx) = val*h;
619  }
620  // y-components
621  for (int j = 0; j < order; j++)
622  for (int i = 0; i <= order; i++)
623  {
624  int idx;
625  if ((idx = dof_map[o++]) < 0)
626  {
627  idx = -1 - idx;
628  }
629 
630  const double h = cp[j+1] - cp[j];
631 
632  double val = 0.0;
633 
634  for (int k = 0; k < nqpt; k++)
635  {
636  const IntegrationPoint &ip1d = ir.IntPoint(k);
637 
638  ip2d.Set2(cp[i], cp[j] + (h*ip1d.x));
639 
640  Trans.SetIntPoint(&ip2d);
641  vc.Eval(xk, Trans, ip2d);
642 
643  // xk^t J tk
644  const double ipval = Trans.Jacobian().InnerProduct(tk + dof2tk[idx]*dim, vk);
645  val += ip1d.weight * ipval;
646  }
647 
648  dofs(idx) = val*h;
649  }
650 }
651 
653  DenseMatrix &shape) const
654 {
655  const int p = order;
656 
657 #ifdef MFEM_THREAD_SAFE
658  Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
659  Vector dshape_cx(p + 1), dshape_cy(p + 1);
660 #endif
661 
663  {
664  cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
665  cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
666  obasis1d.ScaleIntegrated(false);
667  obasis1d.EvalIntegrated(dshape_cx, shape_ox);
668  obasis1d.EvalIntegrated(dshape_cy, shape_oy);
669  }
670  else
671  {
672  cbasis1d.Eval(ip.x, shape_cx);
673  cbasis1d.Eval(ip.y, shape_cy);
674  obasis1d.Eval(ip.x, shape_ox);
675  obasis1d.Eval(ip.y, shape_oy);
676  }
677 
678  int o = 0;
679  // x-components
680  for (int j = 0; j <= p; j++)
681  for (int i = 0; i < p; i++)
682  {
683  int idx, s;
684  if ((idx = dof_map[o++]) < 0)
685  {
686  idx = -1 - idx, s = -1;
687  }
688  else
689  {
690  s = +1;
691  }
692  shape(idx,0) = s*shape_ox(i)*shape_cy(j);
693  shape(idx,1) = 0.;
694  }
695  // y-components
696  for (int j = 0; j < p; j++)
697  for (int i = 0; i <= p; i++)
698  {
699  int idx, s;
700  if ((idx = dof_map[o++]) < 0)
701  {
702  idx = -1 - idx, s = -1;
703  }
704  else
705  {
706  s = +1;
707  }
708  shape(idx,0) = 0.;
709  shape(idx,1) = s*shape_cx(i)*shape_oy(j);
710  }
711 }
712 
714  DenseMatrix &curl_shape) const
715 {
716  const int p = order;
717 
718 #ifdef MFEM_THREAD_SAFE
719  Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
720  Vector dshape_cx(p + 1), dshape_cy(p + 1);
721 #endif
722 
723  cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
724  cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
726  {
727  obasis1d.ScaleIntegrated(false);
728  obasis1d.EvalIntegrated(dshape_cx, shape_ox);
729  obasis1d.EvalIntegrated(dshape_cy, shape_oy);
730  }
731  else
732  {
733  obasis1d.Eval(ip.x, shape_ox);
734  obasis1d.Eval(ip.y, shape_oy);
735  }
736 
737  int o = 0;
738  // x-components
739  for (int j = 0; j <= p; j++)
740  for (int i = 0; i < p; i++)
741  {
742  int idx, s;
743  if ((idx = dof_map[o++]) < 0)
744  {
745  idx = -1 - idx, s = -1;
746  }
747  else
748  {
749  s = +1;
750  }
751  curl_shape(idx,0) = -s*shape_ox(i)*dshape_cy(j);
752  }
753  // y-components
754  for (int j = 0; j < p; j++)
755  for (int i = 0; i <= p; i++)
756  {
757  int idx, s;
758  if ((idx = dof_map[o++]) < 0)
759  {
760  idx = -1 - idx, s = -1;
761  }
762  else
763  {
764  s = +1;
765  }
766  curl_shape(idx,0) = s*dshape_cx(i)*shape_oy(j);
767  }
768 }
769 
770 
771 const double ND_TetrahedronElement::tk[18] =
772 { 1.,0.,0., 0.,1.,0., 0.,0.,1., -1.,1.,0., -1.,0.,1., 0.,-1.,1. };
773 
774 const double ND_TetrahedronElement::c = 1./4.;
775 
777  : VectorFiniteElement(3, Geometry::TETRAHEDRON, p*(p + 2)*(p + 3)/2, p,
778  H_CURL, FunctionSpace::Pk), dof2tk(dof)
779 {
780  const double *eop = poly1d.OpenPoints(p - 1);
781  const double *fop = (p > 1) ? poly1d.OpenPoints(p - 2) : NULL;
782  const double *iop = (p > 2) ? poly1d.OpenPoints(p - 3) : NULL;
783 
784  const int pm1 = p - 1, pm2 = p - 2, pm3 = p - 3;
785 
786 #ifndef MFEM_THREAD_SAFE
787  shape_x.SetSize(p);
788  shape_y.SetSize(p);
789  shape_z.SetSize(p);
790  shape_l.SetSize(p);
791  dshape_x.SetSize(p);
792  dshape_y.SetSize(p);
793  dshape_z.SetSize(p);
794  dshape_l.SetSize(p);
795  u.SetSize(dof, dim);
796 #else
797  Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
798 #endif
799 
800  int o = 0;
801  // edges
802  for (int i = 0; i < p; i++) // (0,1)
803  {
804  Nodes.IntPoint(o).Set3(eop[i], 0., 0.);
805  dof2tk[o++] = 0;
806  }
807  for (int i = 0; i < p; i++) // (0,2)
808  {
809  Nodes.IntPoint(o).Set3(0., eop[i], 0.);
810  dof2tk[o++] = 1;
811  }
812  for (int i = 0; i < p; i++) // (0,3)
813  {
814  Nodes.IntPoint(o).Set3(0., 0., eop[i]);
815  dof2tk[o++] = 2;
816  }
817  for (int i = 0; i < p; i++) // (1,2)
818  {
819  Nodes.IntPoint(o).Set3(eop[pm1-i], eop[i], 0.);
820  dof2tk[o++] = 3;
821  }
822  for (int i = 0; i < p; i++) // (1,3)
823  {
824  Nodes.IntPoint(o).Set3(eop[pm1-i], 0., eop[i]);
825  dof2tk[o++] = 4;
826  }
827  for (int i = 0; i < p; i++) // (2,3)
828  {
829  Nodes.IntPoint(o).Set3(0., eop[pm1-i], eop[i]);
830  dof2tk[o++] = 5;
831  }
832 
833  // faces
834  for (int j = 0; j <= pm2; j++) // (1,2,3)
835  for (int i = 0; i + j <= pm2; i++)
836  {
837  double w = fop[i] + fop[j] + fop[pm2-i-j];
838  Nodes.IntPoint(o).Set3(fop[pm2-i-j]/w, fop[i]/w, fop[j]/w);
839  dof2tk[o++] = 3;
840  Nodes.IntPoint(o).Set3(fop[pm2-i-j]/w, fop[i]/w, fop[j]/w);
841  dof2tk[o++] = 4;
842  }
843  for (int j = 0; j <= pm2; j++) // (0,3,2)
844  for (int i = 0; i + j <= pm2; i++)
845  {
846  double w = fop[i] + fop[j] + fop[pm2-i-j];
847  Nodes.IntPoint(o).Set3(0., fop[j]/w, fop[i]/w);
848  dof2tk[o++] = 2;
849  Nodes.IntPoint(o).Set3(0., fop[j]/w, fop[i]/w);
850  dof2tk[o++] = 1;
851  }
852  for (int j = 0; j <= pm2; j++) // (0,1,3)
853  for (int i = 0; i + j <= pm2; i++)
854  {
855  double w = fop[i] + fop[j] + fop[pm2-i-j];
856  Nodes.IntPoint(o).Set3(fop[i]/w, 0., fop[j]/w);
857  dof2tk[o++] = 0;
858  Nodes.IntPoint(o).Set3(fop[i]/w, 0., fop[j]/w);
859  dof2tk[o++] = 2;
860  }
861  for (int j = 0; j <= pm2; j++) // (0,2,1)
862  for (int i = 0; i + j <= pm2; i++)
863  {
864  double w = fop[i] + fop[j] + fop[pm2-i-j];
865  Nodes.IntPoint(o).Set3(fop[j]/w, fop[i]/w, 0.);
866  dof2tk[o++] = 1;
867  Nodes.IntPoint(o).Set3(fop[j]/w, fop[i]/w, 0.);
868  dof2tk[o++] = 0;
869  }
870 
871  // interior
872  for (int k = 0; k <= pm3; k++)
873  for (int j = 0; j + k <= pm3; j++)
874  for (int i = 0; i + j + k <= pm3; i++)
875  {
876  double w = iop[i] + iop[j] + iop[k] + iop[pm3-i-j-k];
877  Nodes.IntPoint(o).Set3(iop[i]/w, iop[j]/w, iop[k]/w);
878  dof2tk[o++] = 0;
879  Nodes.IntPoint(o).Set3(iop[i]/w, iop[j]/w, iop[k]/w);
880  dof2tk[o++] = 1;
881  Nodes.IntPoint(o).Set3(iop[i]/w, iop[j]/w, iop[k]/w);
882  dof2tk[o++] = 2;
883  }
884 
885  DenseMatrix T(dof);
886  for (int m = 0; m < dof; m++)
887  {
888  const IntegrationPoint &ip = Nodes.IntPoint(m);
889  const double *tm = tk + 3*dof2tk[m];
890  o = 0;
891 
892  poly1d.CalcBasis(pm1, ip.x, shape_x);
893  poly1d.CalcBasis(pm1, ip.y, shape_y);
894  poly1d.CalcBasis(pm1, ip.z, shape_z);
895  poly1d.CalcBasis(pm1, 1. - ip.x - ip.y - ip.z, shape_l);
896 
897  for (int k = 0; k <= pm1; k++)
898  for (int j = 0; j + k <= pm1; j++)
899  for (int i = 0; i + j + k <= pm1; i++)
900  {
901  double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
902  T(o++, m) = s * tm[0];
903  T(o++, m) = s * tm[1];
904  T(o++, m) = s * tm[2];
905  }
906  for (int k = 0; k <= pm1; k++)
907  for (int j = 0; j + k <= pm1; j++)
908  {
909  double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
910  T(o++, m) = s*((ip.y - c)*tm[0] - (ip.x - c)*tm[1]);
911  T(o++, m) = s*((ip.z - c)*tm[0] - (ip.x - c)*tm[2]);
912  }
913  for (int k = 0; k <= pm1; k++)
914  {
915  T(o++, m) =
916  shape_y(pm1-k)*shape_z(k)*((ip.z - c)*tm[1] - (ip.y - c)*tm[2]);
917  }
918  }
919 
920  Ti.Factor(T);
921  // mfem::out << "ND_TetrahedronElement(" << p << ") : "; Ti.TestInversion();
922 }
923 
925  DenseMatrix &shape) const
926 {
927  const int pm1 = order - 1;
928 
929 #ifdef MFEM_THREAD_SAFE
930  const int p = order;
931  Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
932  DenseMatrix u(dof, dim);
933 #endif
934 
935  poly1d.CalcBasis(pm1, ip.x, shape_x);
936  poly1d.CalcBasis(pm1, ip.y, shape_y);
937  poly1d.CalcBasis(pm1, ip.z, shape_z);
938  poly1d.CalcBasis(pm1, 1. - ip.x - ip.y - ip.z, shape_l);
939 
940  int n = 0;
941  for (int k = 0; k <= pm1; k++)
942  for (int j = 0; j + k <= pm1; j++)
943  for (int i = 0; i + j + k <= pm1; i++)
944  {
945  double s = shape_x(i)*shape_y(j)*shape_z(k)*shape_l(pm1-i-j-k);
946  u(n,0) = s; u(n,1) = 0.; u(n,2) = 0.; n++;
947  u(n,0) = 0.; u(n,1) = s; u(n,2) = 0.; n++;
948  u(n,0) = 0.; u(n,1) = 0.; u(n,2) = s; n++;
949  }
950  for (int k = 0; k <= pm1; k++)
951  for (int j = 0; j + k <= pm1; j++)
952  {
953  double s = shape_x(pm1-j-k)*shape_y(j)*shape_z(k);
954  u(n,0) = s*(ip.y - c); u(n,1) = -s*(ip.x - c); u(n,2) = 0.; n++;
955  u(n,0) = s*(ip.z - c); u(n,1) = 0.; u(n,2) = -s*(ip.x - c); n++;
956  }
957  for (int k = 0; k <= pm1; k++)
958  {
959  double s = shape_y(pm1-k)*shape_z(k);
960  u(n,0) = 0.; u(n,1) = s*(ip.z - c); u(n,2) = -s*(ip.y - c); n++;
961  }
962 
963  Ti.Mult(u, shape);
964 }
965 
967  DenseMatrix &curl_shape) const
968 {
969  const int pm1 = order - 1;
970 
971 #ifdef MFEM_THREAD_SAFE
972  const int p = order;
973  Vector shape_x(p), shape_y(p), shape_z(p), shape_l(p);
974  Vector dshape_x(p), dshape_y(p), dshape_z(p), dshape_l(p);
975  DenseMatrix u(dof, dim);
976 #endif
977 
978  poly1d.CalcBasis(pm1, ip.x, shape_x, dshape_x);
979  poly1d.CalcBasis(pm1, ip.y, shape_y, dshape_y);
980  poly1d.CalcBasis(pm1, ip.z, shape_z, dshape_z);
981  poly1d.CalcBasis(pm1, 1. - ip.x - ip.y - ip.z, shape_l, dshape_l);
982 
983  int n = 0;
984  for (int k = 0; k <= pm1; k++)
985  for (int j = 0; j + k <= pm1; j++)
986  for (int i = 0; i + j + k <= pm1; i++)
987  {
988  int l = pm1-i-j-k;
989  const double dx = (dshape_x(i)*shape_l(l) -
990  shape_x(i)*dshape_l(l))*shape_y(j)*shape_z(k);
991  const double dy = (dshape_y(j)*shape_l(l) -
992  shape_y(j)*dshape_l(l))*shape_x(i)*shape_z(k);
993  const double dz = (dshape_z(k)*shape_l(l) -
994  shape_z(k)*dshape_l(l))*shape_x(i)*shape_y(j);
995 
996  u(n,0) = 0.; u(n,1) = dz; u(n,2) = -dy; n++;
997  u(n,0) = -dz; u(n,1) = 0.; u(n,2) = dx; n++;
998  u(n,0) = dy; u(n,1) = -dx; u(n,2) = 0.; n++;
999  }
1000  for (int k = 0; k <= pm1; k++)
1001  for (int j = 0; j + k <= pm1; j++)
1002  {
1003  int i = pm1 - j - k;
1004  // s = shape_x(i)*shape_y(j)*shape_z(k);
1005  // curl of s*(ip.y - c, -(ip.x - c), 0):
1006  u(n,0) = shape_x(i)*(ip.x - c)*shape_y(j)*dshape_z(k);
1007  u(n,1) = shape_x(i)*shape_y(j)*(ip.y - c)*dshape_z(k);
1008  u(n,2) =
1009  -((dshape_x(i)*(ip.x - c) + shape_x(i))*shape_y(j)*shape_z(k) +
1010  (dshape_y(j)*(ip.y - c) + shape_y(j))*shape_x(i)*shape_z(k));
1011  n++;
1012  // curl of s*(ip.z - c, 0, -(ip.x - c)):
1013  u(n,0) = -shape_x(i)*(ip.x - c)*dshape_y(j)*shape_z(k);
1014  u(n,1) = (shape_x(i)*shape_y(j)*(dshape_z(k)*(ip.z - c) + shape_z(k)) +
1015  (dshape_x(i)*(ip.x - c) + shape_x(i))*shape_y(j)*shape_z(k));
1016  u(n,2) = -shape_x(i)*dshape_y(j)*shape_z(k)*(ip.z - c);
1017  n++;
1018  }
1019  for (int k = 0; k <= pm1; k++)
1020  {
1021  int j = pm1 - k;
1022  // curl of shape_y(j)*shape_z(k)*(0, ip.z - c, -(ip.y - c)):
1023  u(n,0) = -((dshape_y(j)*(ip.y - c) + shape_y(j))*shape_z(k) +
1024  shape_y(j)*(dshape_z(k)*(ip.z - c) + shape_z(k)));
1025  u(n,1) = 0.;
1026  u(n,2) = 0.; n++;
1027  }
1028 
1029  Ti.Mult(u, curl_shape);
1030 }
1031 
1032 
1033 const double ND_TriangleElement::tk[8] =
1034 { 1.,0., -1.,1., 0.,-1., 0.,1. };
1035 
1036 const double ND_TriangleElement::c = 1./3.;
1037 
1039  : VectorFiniteElement(2, Geometry::TRIANGLE, p*(p + 2), p,
1040  H_CURL, FunctionSpace::Pk),
1041  dof2tk(dof)
1042 {
1043  const double *eop = poly1d.OpenPoints(p - 1);
1044  const double *iop = (p > 1) ? poly1d.OpenPoints(p - 2) : NULL;
1045 
1046  const int pm1 = p - 1, pm2 = p - 2;
1047 
1048 #ifndef MFEM_THREAD_SAFE
1049  shape_x.SetSize(p);
1050  shape_y.SetSize(p);
1051  shape_l.SetSize(p);
1052  dshape_x.SetSize(p);
1053  dshape_y.SetSize(p);
1054  dshape_l.SetSize(p);
1055  u.SetSize(dof, dim);
1056  curlu.SetSize(dof);
1057 #else
1058  Vector shape_x(p), shape_y(p), shape_l(p);
1059 #endif
1060 
1061  int n = 0;
1062  // edges
1063  for (int i = 0; i < p; i++) // (0,1)
1064  {
1065  Nodes.IntPoint(n).Set2(eop[i], 0.);
1066  dof2tk[n++] = 0;
1067  }
1068  for (int i = 0; i < p; i++) // (1,2)
1069  {
1070  Nodes.IntPoint(n).Set2(eop[pm1-i], eop[i]);
1071  dof2tk[n++] = 1;
1072  }
1073  for (int i = 0; i < p; i++) // (2,0)
1074  {
1075  Nodes.IntPoint(n).Set2(0., eop[pm1-i]);
1076  dof2tk[n++] = 2;
1077  }
1078 
1079  // interior
1080  for (int j = 0; j <= pm2; j++)
1081  for (int i = 0; i + j <= pm2; i++)
1082  {
1083  double w = iop[i] + iop[j] + iop[pm2-i-j];
1084  Nodes.IntPoint(n).Set2(iop[i]/w, iop[j]/w);
1085  dof2tk[n++] = 0;
1086  Nodes.IntPoint(n).Set2(iop[i]/w, iop[j]/w);
1087  dof2tk[n++] = 3;
1088  }
1089 
1090  DenseMatrix T(dof);
1091  for (int m = 0; m < dof; m++)
1092  {
1093  const IntegrationPoint &ip = Nodes.IntPoint(m);
1094  const double *tm = tk + 2*dof2tk[m];
1095  n = 0;
1096 
1097  poly1d.CalcBasis(pm1, ip.x, shape_x);
1098  poly1d.CalcBasis(pm1, ip.y, shape_y);
1099  poly1d.CalcBasis(pm1, 1. - ip.x - ip.y, shape_l);
1100 
1101  for (int j = 0; j <= pm1; j++)
1102  for (int i = 0; i + j <= pm1; i++)
1103  {
1104  double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
1105  T(n++, m) = s * tm[0];
1106  T(n++, m) = s * tm[1];
1107  }
1108  for (int j = 0; j <= pm1; j++)
1109  {
1110  T(n++, m) =
1111  shape_x(pm1-j)*shape_y(j)*((ip.y - c)*tm[0] - (ip.x - c)*tm[1]);
1112  }
1113  }
1114 
1115  Ti.Factor(T);
1116  // mfem::out << "ND_TriangleElement(" << p << ") : "; Ti.TestInversion();
1117 }
1118 
1120  DenseMatrix &shape) const
1121 {
1122  const int pm1 = order - 1;
1123 
1124 #ifdef MFEM_THREAD_SAFE
1125  const int p = order;
1126  Vector shape_x(p), shape_y(p), shape_l(p);
1127  DenseMatrix u(dof, dim);
1128 #endif
1129 
1130  poly1d.CalcBasis(pm1, ip.x, shape_x);
1131  poly1d.CalcBasis(pm1, ip.y, shape_y);
1132  poly1d.CalcBasis(pm1, 1. - ip.x - ip.y, shape_l);
1133 
1134  int n = 0;
1135  for (int j = 0; j <= pm1; j++)
1136  for (int i = 0; i + j <= pm1; i++)
1137  {
1138  double s = shape_x(i)*shape_y(j)*shape_l(pm1-i-j);
1139  u(n,0) = s; u(n,1) = 0; n++;
1140  u(n,0) = 0; u(n,1) = s; n++;
1141  }
1142  for (int j = 0; j <= pm1; j++)
1143  {
1144  double s = shape_x(pm1-j)*shape_y(j);
1145  u(n,0) = s*(ip.y - c);
1146  u(n,1) = -s*(ip.x - c);
1147  n++;
1148  }
1149 
1150  Ti.Mult(u, shape);
1151 }
1152 
1154  DenseMatrix &curl_shape) const
1155 {
1156  const int pm1 = order - 1;
1157 
1158 #ifdef MFEM_THREAD_SAFE
1159  const int p = order;
1160  Vector shape_x(p), shape_y(p), shape_l(p);
1161  Vector dshape_x(p), dshape_y(p), dshape_l(p);
1162  Vector curlu(dof);
1163 #endif
1164 
1165  poly1d.CalcBasis(pm1, ip.x, shape_x, dshape_x);
1166  poly1d.CalcBasis(pm1, ip.y, shape_y, dshape_y);
1167  poly1d.CalcBasis(pm1, 1. - ip.x - ip.y, shape_l, dshape_l);
1168 
1169  int n = 0;
1170  for (int j = 0; j <= pm1; j++)
1171  for (int i = 0; i + j <= pm1; i++)
1172  {
1173  int l = pm1-i-j;
1174  const double dx = (dshape_x(i)*shape_l(l) -
1175  shape_x(i)*dshape_l(l)) * shape_y(j);
1176  const double dy = (dshape_y(j)*shape_l(l) -
1177  shape_y(j)*dshape_l(l)) * shape_x(i);
1178 
1179  curlu(n++) = -dy;
1180  curlu(n++) = dx;
1181  }
1182 
1183  for (int j = 0; j <= pm1; j++)
1184  {
1185  int i = pm1 - j;
1186  // curl of shape_x(i)*shape_y(j) * (ip.y - c, -(ip.x - c), 0):
1187  curlu(n++) = -((dshape_x(i)*(ip.x - c) + shape_x(i)) * shape_y(j) +
1188  (dshape_y(j)*(ip.y - c) + shape_y(j)) * shape_x(i));
1189  }
1190 
1191  Vector curl2d(curl_shape.Data(),dof);
1192  Ti.Mult(curlu, curl2d);
1193 }
1194 
1195 
1196 const double ND_SegmentElement::tk[1] = { 1. };
1197 
1198 ND_SegmentElement::ND_SegmentElement(const int p, const int ob_type)
1199  : VectorTensorFiniteElement(1, p, p - 1, ob_type, H_CURL,
1200  DofMapType::L2_DOF_MAP),
1201  dof2tk(dof)
1202 {
1203  if (obasis1d.IsIntegratedType()) { is_nodal = false; }
1204 
1205  const double *op = poly1d.OpenPoints(p - 1, ob_type);
1206 
1207  // set dof2tk and Nodes
1208  for (int i = 0; i < p; i++)
1209  {
1210  dof2tk[i] = 0;
1211  Nodes.IntPoint(i).x = op[i];
1212  }
1213 }
1214 
1216  DenseMatrix &shape) const
1217 {
1218  Vector vshape(shape.Data(), dof);
1219 
1220  obasis1d.Eval(ip.x, vshape);
1221 }
1222 
1223 const double ND_WedgeElement::tk[15] =
1224 { 1.,0.,0., -1.,1.,0., 0.,-1.,0., 0.,0.,1., 0.,1.,0. };
1225 
1227  const int cb_type,
1228  const int ob_type)
1229  : VectorFiniteElement(3, Geometry::PRISM,
1230  3 * p * ((p + 1) * (p + 2))/2, p,
1231  H_CURL, FunctionSpace::Qk),
1232  dof2tk(dof),
1233  t_dof(dof),
1234  s_dof(dof),
1235  H1TriangleFE(p, cb_type),
1236  NDTriangleFE(p),
1237  H1SegmentFE(p, cb_type),
1238  NDSegmentFE(p, ob_type)
1239 {
1240  MFEM_ASSERT(H1TriangleFE.GetDof() * NDSegmentFE.GetDof() +
1241  NDTriangleFE.GetDof() * H1SegmentFE.GetDof() == dof,
1242  "Mismatch in number of degrees of freedom "
1243  "when building ND_WedgeElement!");
1244 
1245 #ifndef MFEM_THREAD_SAFE
1246  t1_shape.SetSize(H1TriangleFE.GetDof());
1247  s1_shape.SetSize(H1SegmentFE.GetDof());
1248  tn_shape.SetSize(NDTriangleFE.GetDof(), 2);
1249  sn_shape.SetSize(NDSegmentFE.GetDof(), 1);
1250  t1_dshape.SetSize(H1TriangleFE.GetDof(), 2);
1251  s1_dshape.SetSize(H1SegmentFE.GetDof(), 1);
1252  tn_dshape.SetSize(NDTriangleFE.GetDof(), 1);
1253 #endif
1254 
1255  const int pm1 = p - 1, pm2 = p - 2;
1256 
1257  const IntegrationRule &t1_n = H1TriangleFE.GetNodes();
1258  const IntegrationRule &tn_n = NDTriangleFE.GetNodes();
1259  const IntegrationRule &s1_n = H1SegmentFE.GetNodes();
1260  const IntegrationRule &sn_n = NDSegmentFE.GetNodes();
1261 
1262  // edges
1263  int o = 0;
1264  for (int i = 0; i < p; i++) // (0,1)
1265  {
1266  t_dof[o] = i; s_dof[o] = 0; dof2tk[o] = 0;
1267  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1268  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1269  o++;
1270  }
1271  for (int i = 0; i < p; i++) // (1,2)
1272  {
1273  t_dof[o] = p + i; s_dof[o] = 0; dof2tk[o] = 1;
1274  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1275  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1276  o++;
1277  }
1278  for (int i = 0; i < p; i++) // (2,0)
1279  {
1280  t_dof[o] = 2 * p + i; s_dof[o] = 0; dof2tk[o] = 2;
1281  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1282  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1283  o++;
1284  }
1285  for (int i = 0; i < p; i++) // (3,4)
1286  {
1287  t_dof[o] = i; s_dof[o] = 1; dof2tk[o] = 0;
1288  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1289  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1290  o++;
1291  }
1292  for (int i = 0; i < p; i++) // (4,5)
1293  {
1294  t_dof[o] = p + i; s_dof[o] = 1; dof2tk[o] = 1;
1295  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1296  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1297  o++;
1298  }
1299  for (int i = 0; i < p; i++) // (5,3)
1300  {
1301  t_dof[o] = 2 * p + i; s_dof[o] = 1; dof2tk[o] = 2;
1302  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1303  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1304  o++;
1305  }
1306  for (int i = 0; i < p; i++) // (0,3)
1307  {
1308  t_dof[o] = 0; s_dof[o] = i; dof2tk[o] = 3;
1309  const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1310  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1311  o++;
1312  }
1313  for (int i = 0; i < p; i++) // (1,4)
1314  {
1315  t_dof[o] = 1; s_dof[o] = i; dof2tk[o] = 3;
1316  const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1317  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1318  o++;
1319  }
1320  for (int i = 0; i < p; i++) // (2,5)
1321  {
1322  t_dof[o] = 2; s_dof[o] = i; dof2tk[o] = 3;
1323  const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1324  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1325  o++;
1326  }
1327 
1328  // faces
1329  // (0,2,1) -- bottom
1330  int l = 0;
1331  for (int j = 0; j <= pm2; j++)
1332  for (int i = 0; i + j <= pm2; i++)
1333  {
1334  l = j + ( 2 * p - 1 - i) * i / 2;
1335  t_dof[o] = 3 * p + 2*l+1; s_dof[o] = 0; dof2tk[o] = 4;
1336  const IntegrationPoint & t_ip0 = tn_n.IntPoint(t_dof[o]);
1337  Nodes.IntPoint(o).Set3(t_ip0.x, t_ip0.y, s1_n.IntPoint(s_dof[o]).x);
1338  o++;
1339  t_dof[o] = 3 * p + 2*l; s_dof[o] = 0; dof2tk[o] = 0;
1340  const IntegrationPoint & t_ip1 = tn_n.IntPoint(t_dof[o]);
1341  Nodes.IntPoint(o).Set3(t_ip1.x, t_ip1.y, s1_n.IntPoint(s_dof[o]).x);
1342  o++;
1343  }
1344  // (3,4,5) -- top
1345  int m = 0;
1346  for (int j = 0; j <= pm2; j++)
1347  for (int i = 0; i + j <= pm2; i++)
1348  {
1349  t_dof[o] = 3 * p + m; s_dof[o] = 1; dof2tk[o] = 0; m++;
1350  const IntegrationPoint & t_ip0 = tn_n.IntPoint(t_dof[o]);
1351  Nodes.IntPoint(o).Set3(t_ip0.x, t_ip0.y, s1_n.IntPoint(s_dof[o]).x);
1352  o++;
1353  t_dof[o] = 3 * p + m; s_dof[o] = 1; dof2tk[o] = 4; m++;
1354  const IntegrationPoint & t_ip1 = tn_n.IntPoint(t_dof[o]);
1355  Nodes.IntPoint(o).Set3(t_ip1.x, t_ip1.y, s1_n.IntPoint(s_dof[o]).x);
1356  o++;
1357  }
1358  // (0, 1, 4, 3) -- xz plane
1359  for (int j = 2; j <= p; j++)
1360  for (int i = 0; i < p; i++)
1361  {
1362  t_dof[o] = i; s_dof[o] = j; dof2tk[o] = 0;
1363  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1364  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1365  o++;
1366  }
1367  for (int j = 0; j < p; j++)
1368  for (int i = 0; i < pm1; i++)
1369  {
1370  t_dof[o] = 3 + i; s_dof[o] = j; dof2tk[o] = 3;
1371  const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1372  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1373  o++;
1374  }
1375  // (1, 2, 5, 4) -- (y-x)z plane
1376  for (int j = 2; j <= p; j++)
1377  for (int i = 0; i < p; i++)
1378  {
1379  t_dof[o] = p + i; s_dof[o] = j; dof2tk[o] = 1;
1380  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1381  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1382  o++;
1383  }
1384  for (int j = 0; j < p; j++)
1385  for (int i = 0; i < pm1; i++)
1386  {
1387  t_dof[o] = p + 2 + i; s_dof[o] = j; dof2tk[o] = 3;
1388  const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1389  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1390  o++;
1391  }
1392  // (2, 0, 3, 5) -- yz plane
1393  for (int j = 2; j <= p; j++)
1394  for (int i = 0; i < p; i++)
1395  {
1396  t_dof[o] = 2 * p + i; s_dof[o] = j; dof2tk[o] = 2;
1397  const IntegrationPoint & t_ip = tn_n.IntPoint(t_dof[o]);
1398  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, s1_n.IntPoint(s_dof[o]).x);
1399  o++;
1400  }
1401  for (int j = 0; j < p; j++)
1402  for (int i = 0; i < pm1; i++)
1403  {
1404  t_dof[o] = 2 * p + 1 + i; s_dof[o] = j; dof2tk[o] = 3;
1405  const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1406  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1407  o++;
1408  }
1409 
1410  // interior
1411  for (int k = 2; k <= p; k++)
1412  {
1413  l = 0;
1414  for (int j = 0; j <= pm2; j++)
1415  for (int i = 0; i + j <= pm2; i++)
1416  {
1417  t_dof[o] = 3 * p + l; s_dof[o] = k; dof2tk[o] = 0; l++;
1418  const IntegrationPoint & t_ip0 = tn_n.IntPoint(t_dof[o]);
1419  Nodes.IntPoint(o).Set3(t_ip0.x, t_ip0.y, s1_n.IntPoint(s_dof[o]).x);
1420  o++;
1421  t_dof[o] = 3 * p + l; s_dof[o] = k; dof2tk[o] = 4; l++;
1422  const IntegrationPoint & t_ip1 = tn_n.IntPoint(t_dof[o]);
1423  Nodes.IntPoint(o).Set3(t_ip1.x, t_ip1.y, s1_n.IntPoint(s_dof[o]).x);
1424  o++;
1425  }
1426  }
1427  for (int k = 0; k < p; k++)
1428  {
1429  l = 0;
1430  for (int j = 0; j < pm2; j++)
1431  for (int i = 0; i + j < pm2; i++)
1432  {
1433  t_dof[o] = 3 * p + l; s_dof[o] = k; dof2tk[o] = 3; l++;
1434  const IntegrationPoint & t_ip = t1_n.IntPoint(t_dof[o]);
1435  Nodes.IntPoint(o).Set3(t_ip.x, t_ip.y, sn_n.IntPoint(s_dof[o]).x);
1436  o++;
1437  }
1438  }
1439 }
1440 
1442  DenseMatrix &shape) const
1443 {
1444 #ifdef MFEM_THREAD_SAFE
1445  Vector t1_shape(H1TriangleFE.GetDof());
1446  Vector s1_shape(H1SegmentFE.GetDof());
1447  DenseMatrix tn_shape(NDTriangleFE.GetDof(), 2);
1448  DenseMatrix sn_shape(NDSegmentFE.GetDof(), 1);
1449 #endif
1450 
1451  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1452 
1453  H1TriangleFE.CalcShape(ip, t1_shape);
1454  NDTriangleFE.CalcVShape(ip, tn_shape);
1455  H1SegmentFE.CalcShape(ipz, s1_shape);
1456  NDSegmentFE.CalcVShape(ipz, sn_shape);
1457 
1458  for (int i=0; i<dof; i++)
1459  {
1460  if ( dof2tk[i] != 3 )
1461  {
1462  shape(i, 0) = tn_shape(t_dof[i], 0) * s1_shape[s_dof[i]];
1463  shape(i, 1) = tn_shape(t_dof[i], 1) * s1_shape[s_dof[i]];
1464  shape(i, 2) = 0.0;
1465  }
1466  else
1467  {
1468  shape(i, 0) = 0.0;
1469  shape(i, 1) = 0.0;
1470  shape(i, 2) = t1_shape[t_dof[i]] * sn_shape(s_dof[i], 0);
1471  }
1472  }
1473 }
1474 
1476  DenseMatrix &curl_shape) const
1477 {
1478 #ifdef MFEM_THREAD_SAFE
1479  Vector s1_shape(H1SegmentFE.GetDof());
1480  DenseMatrix t1_dshape(H1TriangleFE.GetDof(), 2);
1481  DenseMatrix s1_dshape(H1SegmentFE.GetDof(), 1);
1482  DenseMatrix tn_shape(NDTriangleFE.GetDof(), 2);
1483  DenseMatrix sn_shape(NDSegmentFE.GetDof(), 1);
1484  DenseMatrix tn_dshape(NDTriangleFE.GetDof(), 1);
1485 #endif
1486 
1487  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1488 
1489  H1TriangleFE.CalcDShape(ip, t1_dshape);
1490  H1SegmentFE.CalcShape(ipz, s1_shape);
1491  H1SegmentFE.CalcDShape(ipz, s1_dshape);
1492  NDTriangleFE.CalcVShape(ip, tn_shape);
1493  NDTriangleFE.CalcCurlShape(ip, tn_dshape);
1494  NDSegmentFE.CalcVShape(ipz, sn_shape);
1495 
1496  for (int i=0; i<dof; i++)
1497  {
1498  if ( dof2tk[i] != 3 )
1499  {
1500  curl_shape(i, 0) = -tn_shape(t_dof[i], 1) * s1_dshape(s_dof[i], 0);
1501  curl_shape(i, 1) = tn_shape(t_dof[i], 0) * s1_dshape(s_dof[i], 0);
1502  curl_shape(i, 2) = tn_dshape(t_dof[i], 0) * s1_shape[s_dof[i]];
1503  }
1504  else
1505  {
1506  curl_shape(i, 0) = t1_dshape(t_dof[i], 1) * sn_shape(s_dof[i], 0);
1507  curl_shape(i, 1) = -t1_dshape(t_dof[i], 0) * sn_shape(s_dof[i], 0);
1508  curl_shape(i, 2) = 0.0;
1509  }
1510  }
1511 }
1512 
1514  : VectorFiniteElement(1, Geometry::POINT, 2, p,
1515  H_CURL, FunctionSpace::Pk)
1516 {
1517  // VectorFiniteElement::SetDerivMembers doesn't support 0D H_CURL elements
1518  // so we mimic a 1D element and then correct the dimension here.
1519  dim = 0;
1520  vdim = 2;
1521  cdim = 0;
1522 }
1523 
1525  DenseMatrix &shape) const
1526 {
1527  shape(0,0) = 1.0;
1528  shape(0,1) = 0.0;
1529 
1530  shape(1,0) = 0.0;
1531  shape(1,1) = 1.0;
1532 }
1533 
1535  DenseMatrix &shape) const
1536 {
1537  CalcVShape(Trans.GetIntPoint(), shape);
1538 }
1539 
1540 const double ND_R1D_SegmentElement::tk[9] = { 1.,0.,0., 0.,1.,0., 0.,0.,1. };
1541 
1543  const int cb_type,
1544  const int ob_type)
1545  : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 2, p,
1546  H_CURL, FunctionSpace::Pk),
1547  dof2tk(dof),
1548  cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
1549  obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
1550 {
1551  // Override default types for VectorFiniteElements
1552  deriv_type = CURL;
1555 
1556  // Override default dimensions for VectorFiniteElements
1557  vdim = 3;
1558  cdim = 3;
1559 
1560  const double *cp = poly1d.ClosedPoints(p, cb_type);
1561  const double *op = poly1d.OpenPoints(p - 1, ob_type);
1562 
1563 #ifndef MFEM_THREAD_SAFE
1564  shape_cx.SetSize(p + 1);
1565  shape_ox.SetSize(p);
1566  dshape_cx.SetSize(p + 1);
1567 #endif
1568 
1569  dof_map.SetSize(dof);
1570 
1571  int o = 0;
1572  // nodes
1573  // (0)
1574  Nodes.IntPoint(o).x = cp[0]; // y-directed
1575  dof_map[p] = o; dof2tk[o++] = 1;
1576  Nodes.IntPoint(o).x = cp[0]; // z-directed
1577  dof_map[2*p+1] = o; dof2tk[o++] = 2;
1578 
1579  // (1)
1580  Nodes.IntPoint(o).x = cp[p]; // y-directed
1581  dof_map[2*p] = o; dof2tk[o++] = 1;
1582  Nodes.IntPoint(o).x = cp[p]; // z-directed
1583  dof_map[3*p+1] = o; dof2tk[o++] = 2;
1584 
1585  // interior
1586  // x-components
1587  for (int i = 0; i < p; i++)
1588  {
1589  Nodes.IntPoint(o).x = op[i];
1590  dof_map[i] = o; dof2tk[o++] = 0;
1591  }
1592  // y-components
1593  for (int i = 1; i < p; i++)
1594  {
1595  Nodes.IntPoint(o).x = cp[i];
1596  dof_map[p+i] = o; dof2tk[o++] = 1;
1597  }
1598  // z-components
1599  for (int i = 1; i < p; i++)
1600  {
1601  Nodes.IntPoint(o).x = cp[i];
1602  dof_map[2*p+1+i] = o; dof2tk[o++] = 2;
1603  }
1604 }
1605 
1607  DenseMatrix &shape) const
1608 {
1609  const int p = order;
1610 
1611 #ifdef MFEM_THREAD_SAFE
1612  Vector shape_cx(p + 1), shape_ox(p);
1613 #endif
1614 
1615  cbasis1d.Eval(ip.x, shape_cx);
1616  obasis1d.Eval(ip.x, shape_ox);
1617 
1618  int o = 0;
1619  // x-components
1620  for (int i = 0; i < p; i++)
1621  {
1622  int idx = dof_map[o++];
1623  shape(idx,0) = shape_ox(i);
1624  shape(idx,1) = 0.;
1625  shape(idx,2) = 0.;
1626  }
1627  // y-components
1628  for (int i = 0; i <= p; i++)
1629  {
1630  int idx = dof_map[o++];
1631  shape(idx,0) = 0.;
1632  shape(idx,1) = shape_cx(i);
1633  shape(idx,2) = 0.;
1634  }
1635  // z-components
1636  for (int i = 0; i <= p; i++)
1637  {
1638  int idx = dof_map[o++];
1639  shape(idx,0) = 0.;
1640  shape(idx,1) = 0.;
1641  shape(idx,2) = shape_cx(i);
1642  }
1643 }
1644 
1646  DenseMatrix &shape) const
1647 {
1648  CalcVShape(Trans.GetIntPoint(), shape);
1649  const DenseMatrix & JI = Trans.InverseJacobian();
1650  MFEM_ASSERT(JI.Width() == 1 && JI.Height() == 1,
1651  "ND_R1D_SegmentElement cannot be embedded in "
1652  "2 or 3 dimensional spaces");
1653  for (int i=0; i<dof; i++)
1654  {
1655  shape(i, 0) *= JI(0,0);
1656  }
1657 }
1658 
1660  DenseMatrix &curl_shape) const
1661 {
1662  const int p = order;
1663 
1664 #ifdef MFEM_THREAD_SAFE
1665  Vector shape_cx(p + 1), shape_ox(p);
1666  Vector dshape_cx(p + 1);
1667 #endif
1668 
1669  cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
1670  obasis1d.Eval(ip.x, shape_ox);
1671 
1672  int o = 0;
1673  // x-components
1674  for (int i = 0; i < p; i++)
1675  {
1676  int idx = dof_map[o++];
1677  curl_shape(idx,0) = 0.;
1678  curl_shape(idx,1) = 0.;
1679  curl_shape(idx,2) = 0.;
1680  }
1681  // y-components
1682  for (int i = 0; i <= p; i++)
1683  {
1684  int idx = dof_map[o++];
1685  curl_shape(idx,0) = 0.;
1686  curl_shape(idx,1) = 0.;
1687  curl_shape(idx,2) = dshape_cx(i);
1688  }
1689  // z-components
1690  for (int i = 0; i <= p; i++)
1691  {
1692  int idx = dof_map[o++];
1693  curl_shape(idx,0) = 0.;
1694  curl_shape(idx,1) = -dshape_cx(i);
1695  curl_shape(idx,2) = 0.;
1696  }
1697 }
1698 
1700  DenseMatrix &curl_shape) const
1701 {
1702  CalcCurlShape(Trans.GetIntPoint(), curl_shape);
1703  const DenseMatrix & J = Trans.Jacobian();
1704  MFEM_ASSERT(J.Width() == 1 && J.Height() == 1,
1705  "ND_R1D_SegmentElement cannot be embedded in "
1706  "2 or 3 dimensional spaces");
1707  curl_shape *= (1.0 / J.Weight());
1708 }
1709 
1712  Vector &dofs) const
1713 {
1714  double data[3];
1715  Vector vk(data, 3);
1716 
1717  for (int k = 0; k < dof; k++)
1718  {
1719  Trans.SetIntPoint(&Nodes.IntPoint(k));
1720 
1721  vc.Eval(vk, Trans, Nodes.IntPoint(k));
1722  // dof_k = vk^t J tk
1723  Vector t(const_cast<double*>(&tk[dof2tk[k] * 3]), 3);
1724  dofs(k) = Trans.Jacobian()(0,0) * t(0) * vk(0) +
1725  t(1) * vk(1) + t(2) * vk(2);
1726  }
1727 
1728 }
1729 
1732  DenseMatrix &I) const
1733 {
1734  if (fe.GetRangeType() == SCALAR)
1735  {
1736  double vk[Geometry::MaxDim];
1737  Vector shape(fe.GetDof());
1738 
1739  double * tk_ptr = const_cast<double*>(tk);
1740 
1741  I.SetSize(dof, vdim*fe.GetDof());
1742  for (int k = 0; k < dof; k++)
1743  {
1744  const IntegrationPoint &ip = Nodes.IntPoint(k);
1745 
1746  Vector t1(&tk_ptr[dof2tk[k] * 3], 1);
1747  Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
1748 
1749  fe.CalcShape(ip, shape);
1750  Trans.SetIntPoint(&ip);
1751  // Transform ND edge tengents from reference to physical space
1752  // vk = J tk
1753  Trans.Jacobian().Mult(t1, vk);
1754  vk[1] = t3[1];
1755  vk[2] = t3[2];
1756  if (fe.GetMapType() == INTEGRAL)
1757  {
1758  double w = 1.0/Trans.Weight();
1759  for (int d = 0; d < vdim; d++)
1760  {
1761  vk[d] *= w;
1762  }
1763  }
1764 
1765  for (int j = 0; j < shape.Size(); j++)
1766  {
1767  double s = shape(j);
1768  if (fabs(s) < 1e-12)
1769  {
1770  s = 0.0;
1771  }
1772  // Project scalar basis function multiplied by each coordinate
1773  // direction onto the transformed edge tangents
1774  for (int d = 0; d < vdim; d++)
1775  {
1776  I(k, j + d*shape.Size()) = s*vk[d];
1777  }
1778  }
1779  }
1780  }
1781  else
1782  {
1783  double vk[Geometry::MaxDim];
1784  DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
1785 
1786  double * tk_ptr = const_cast<double*>(tk);
1787 
1788  I.SetSize(dof, fe.GetDof());
1789  for (int k = 0; k < dof; k++)
1790  {
1791  const IntegrationPoint &ip = Nodes.IntPoint(k);
1792 
1793  Vector t1(&tk_ptr[dof2tk[k] * 3], 1);
1794  Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
1795 
1796  Trans.SetIntPoint(&ip);
1797  // Transform ND edge tangents from reference to physical space
1798  // vk = J tk
1799  Trans.Jacobian().Mult(t1, vk);
1800  // Compute fe basis functions in physical space
1801  fe.CalcVShape(Trans, vshape);
1802  // Project fe basis functions onto transformed edge tangents
1803  for (int j=0; j<vshape.Height(); j++)
1804  {
1805  I(k, j) = 0.0;
1806  I(k, j) += vshape(j, 0) * vk[0];
1807  if (vshape.Width() == 3)
1808  {
1809  I(k, j) += vshape(j, 1) * t3(1);
1810  I(k, j) += vshape(j, 2) * t3(2);
1811  }
1812  }
1813  }
1814  }
1815 }
1816 
1817 const double ND_R2D_SegmentElement::tk[4] = { 1.,0., 0.,1. };
1818 
1820  const int cb_type,
1821  const int ob_type)
1822  : VectorFiniteElement(1, Geometry::SEGMENT, 2 * p + 1, p,
1823  H_CURL, FunctionSpace::Pk),
1824  dof2tk(dof),
1825  cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
1826  obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
1827 {
1828  // Override default dimensions for VectorFiniteElements
1829  vdim = 2;
1830  cdim = 1;
1831 
1832  const double *cp = poly1d.ClosedPoints(p, cb_type);
1833  const double *op = poly1d.OpenPoints(p - 1, ob_type);
1834 
1835 #ifndef MFEM_THREAD_SAFE
1836  shape_cx.SetSize(p + 1);
1837  shape_ox.SetSize(p);
1838  dshape_cx.SetSize(p + 1);
1839 #endif
1840 
1841  dof_map.SetSize(dof);
1842 
1843  int o = 0;
1844  // nodes
1845  // (0)
1846  Nodes.IntPoint(o).x = cp[0]; // z-directed
1847  dof_map[p] = o; dof2tk[o++] = 1;
1848 
1849  // (1)
1850  Nodes.IntPoint(o).x = cp[p]; // z-directed
1851  dof_map[2*p] = o; dof2tk[o++] = 1;
1852 
1853  // interior
1854  // x-components
1855  for (int i = 0; i < p; i++)
1856  {
1857  Nodes.IntPoint(o).x = op[i];
1858  dof_map[i] = o; dof2tk[o++] = 0;
1859  }
1860  // z-components
1861  for (int i = 1; i < p; i++)
1862  {
1863  Nodes.IntPoint(o).x = cp[i];
1864  dof_map[p+i] = o; dof2tk[o++] = 1;
1865  }
1866 }
1867 
1869  DenseMatrix &shape) const
1870 {
1871  const int p = order;
1872 
1873 #ifdef MFEM_THREAD_SAFE
1874  Vector shape_cx(p + 1), shape_ox(p);
1875 #endif
1876 
1877  cbasis1d.Eval(ip.x, shape_cx);
1878  obasis1d.Eval(ip.x, shape_ox);
1879 
1880  int o = 0;
1881  // x-components
1882  for (int i = 0; i < p; i++)
1883  {
1884  int idx = dof_map[o++];
1885  shape(idx,0) = shape_ox(i);
1886  shape(idx,1) = 0.;
1887  }
1888  // z-components
1889  for (int i = 0; i <= p; i++)
1890  {
1891  int idx = dof_map[o++];
1892  shape(idx,0) = 0.;
1893  shape(idx,1) = shape_cx(i);
1894  }
1895 }
1896 
1898  DenseMatrix &shape) const
1899 {
1900  CalcVShape(Trans.GetIntPoint(), shape);
1901  const DenseMatrix & JI = Trans.InverseJacobian();
1902  MFEM_ASSERT(JI.Width() == 1 && JI.Height() == 1,
1903  "ND_R2D_SegmentElement cannot be embedded in "
1904  "2 or 3 dimensional spaces");
1905  for (int i=0; i<dof; i++)
1906  {
1907  shape(i, 0) *= JI(0,0);
1908  }
1909 }
1910 
1912  DenseMatrix &curl_shape) const
1913 {
1914  const int p = order;
1915 
1916 #ifdef MFEM_THREAD_SAFE
1917  Vector shape_cx(p + 1), shape_ox(p);
1918  Vector dshape_cx(p + 1);
1919 #endif
1920 
1921  cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
1922  obasis1d.Eval(ip.x, shape_ox);
1923 
1924  int o = 0;
1925  // x-components
1926  for (int i = 0; i < p; i++)
1927  {
1928  int idx = dof_map[o++];
1929  curl_shape(idx,0) = 0.;
1930  }
1931  // z-components
1932  for (int i = 0; i <= p; i++)
1933  {
1934  int idx = dof_map[o++];
1935  curl_shape(idx,0) = -dshape_cx(i);
1936  }
1937 }
1938 
1939 void ND_R2D_SegmentElement::LocalInterpolation(const VectorFiniteElement &cfe,
1941  DenseMatrix &I) const
1942 {
1943  double vk[Geometry::MaxDim]; vk[1] = 0.0; vk[2] = 0.0;
1944  Vector xk(vk, dim);
1945  IntegrationPoint ip;
1946  DenseMatrix vshape(cfe.GetDof(), vdim);
1947 
1948  double * tk_ptr = const_cast<double*>(tk);
1949 
1950  I.SetSize(dof, vshape.Height());
1951 
1952  // assuming Trans is linear; this should be ok for all refinement types
1954  const DenseMatrix &J = Trans.Jacobian();
1955  for (int k = 0; k < dof; k++)
1956  {
1957  Vector t1(&tk_ptr[dof2tk[k] * 2], 1);
1958  Vector t2(&tk_ptr[dof2tk[k] * 2], 2);
1959 
1960  Trans.Transform(Nodes.IntPoint(k), xk);
1961  ip.Set3(vk);
1962  cfe.CalcVShape(ip, vshape);
1963  // xk = J t_k
1964  J.Mult(t1, vk);
1965  // I_k = vshape_k.J.t_k, k=1,...,Dof
1966  for (int j = 0; j < vshape.Height(); j++)
1967  {
1968  double Ikj = 0.;
1969  for (int i = 0; i < dim; i++)
1970  {
1971  Ikj += vshape(j, i) * vk[i];
1972  }
1973  Ikj += vshape(j, 1) * t2(1);
1974  I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1975  }
1976  }
1977 }
1978 
1980  ElementTransformation &Trans,
1981  Vector &dofs) const
1982 {
1983  double data[3];
1984  Vector vk1(data, 1);
1985  Vector vk2(data, 2);
1986  Vector vk3(data, 3);
1987 
1988  double * tk_ptr = const_cast<double*>(tk);
1989 
1990  for (int k = 0; k < dof; k++)
1991  {
1992  Trans.SetIntPoint(&Nodes.IntPoint(k));
1993 
1994  vc.Eval(vk3, Trans, Nodes.IntPoint(k));
1995  // dof_k = vk^t J tk
1996  Vector t1(&tk_ptr[dof2tk[k] * 2], 1);
1997  Vector t2(&tk_ptr[dof2tk[k] * 2], 2);
1998 
1999  dofs(k) = Trans.Jacobian().InnerProduct(t1, vk2) + t2(1) * vk3(2);
2000  }
2001 
2002 }
2003 
2005  const double *tk_fe)
2006  : VectorFiniteElement(2, G, Do, p,
2007  H_CURL, FunctionSpace::Pk),
2008  tk(tk_fe),
2009  dof_map(dof),
2010  dof2tk(dof)
2011 {
2012  // Override default types for VectorFiniteElements
2013  deriv_type = CURL;
2016 
2017  // Override default dimensions for VectorFiniteElements
2018  vdim = 3;
2019  cdim = 3;
2020 }
2021 
2023  DenseMatrix &shape) const
2024 {
2025  CalcVShape(Trans.GetIntPoint(), shape);
2026  const DenseMatrix & JI = Trans.InverseJacobian();
2027  MFEM_ASSERT(JI.Width() == 2 && JI.Height() == 2,
2028  "ND_R2D_FiniteElement cannot be embedded in "
2029  "3 dimensional spaces");
2030  for (int i=0; i<dof; i++)
2031  {
2032  double sx = shape(i, 0);
2033  double sy = shape(i, 1);
2034  shape(i, 0) = sx * JI(0, 0) + sy * JI(1, 0);
2035  shape(i, 1) = sx * JI(0, 1) + sy * JI(1, 1);
2036  }
2037 }
2038 
2040  DenseMatrix &curl_shape) const
2041 {
2042  CalcCurlShape(Trans.GetIntPoint(), curl_shape);
2043  const DenseMatrix & J = Trans.Jacobian();
2044  MFEM_ASSERT(J.Width() == 2 && J.Height() == 2,
2045  "ND_R2D_FiniteElement cannot be embedded in "
2046  "3 dimensional spaces");
2047  for (int i=0; i<dof; i++)
2048  {
2049  double sx = curl_shape(i, 0);
2050  double sy = curl_shape(i, 1);
2051  curl_shape(i, 0) = sx * J(0, 0) + sy * J(0, 1);
2052  curl_shape(i, 1) = sx * J(1, 0) + sy * J(1, 1);
2053  }
2054  curl_shape *= (1.0 / Trans.Weight());
2055 }
2056 
2057 void ND_R2D_FiniteElement::LocalInterpolation(
2058  const VectorFiniteElement &cfe,
2059  ElementTransformation &Trans,
2060  DenseMatrix &I) const
2061 {
2062  double vk[Geometry::MaxDim]; vk[2] = 0.0;
2063  Vector xk(vk, dim);
2064  IntegrationPoint ip;
2065 #ifdef MFEM_THREAD_SAFE
2066  DenseMatrix vshape(cfe.GetDof(), vdim);
2067 #else
2068  vshape.SetSize(cfe.GetDof(), vdim);
2069 #endif
2070 
2071  double * tk_ptr = const_cast<double*>(tk);
2072 
2073  I.SetSize(dof, vshape.Height());
2074 
2075  // assuming Trans is linear; this should be ok for all refinement types
2077  const DenseMatrix &J = Trans.Jacobian();
2078  for (int k = 0; k < dof; k++)
2079  {
2080  Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
2081  Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2082 
2083  Trans.Transform(Nodes.IntPoint(k), xk);
2084  ip.Set3(vk);
2085  cfe.CalcVShape(ip, vshape);
2086  // xk = J t_k
2087  J.Mult(t2, vk);
2088  // I_k = vshape_k.J.t_k, k=1,...,Dof
2089  for (int j = 0; j < vshape.Height(); j++)
2090  {
2091  double Ikj = 0.;
2092  for (int i = 0; i < dim; i++)
2093  {
2094  Ikj += vshape(j, i) * vk[i];
2095  }
2096  Ikj += vshape(j, 2) * t3(2);
2097  I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
2098  }
2099  }
2100 }
2101 
2103  DenseMatrix &R) const
2104 {
2105  double pt_data[Geometry::MaxDim];
2106  IntegrationPoint ip;
2107  Vector pt(pt_data, dim);
2108 
2109 #ifdef MFEM_THREAD_SAFE
2110  DenseMatrix vshape(dof, vdim);
2111 #endif
2112 
2113  double * tk_ptr = const_cast<double*>(tk);
2114 
2116  const DenseMatrix &Jinv = Trans.InverseJacobian();
2117  for (int j = 0; j < dof; j++)
2118  {
2119  Vector t2(&tk_ptr[dof2tk[j] * 3], 2);
2120  Vector t3(&tk_ptr[dof2tk[j] * 3], 3);
2121 
2122  InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
2123  ip.Set(pt_data, dim);
2124  if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
2125  {
2126  CalcVShape(ip, vshape);
2127  Jinv.Mult(t2, pt_data);
2128  for (int k = 0; k < dof; k++)
2129  {
2130  double R_jk = 0.0;
2131  for (int d = 0; d < dim; d++)
2132  {
2133  R_jk += vshape(k,d)*pt_data[d];
2134  }
2135  R_jk += vshape(k, 2) * t3(2);
2136  R(j,k) = R_jk;
2137  }
2138  }
2139  else
2140  {
2141  // Set the whole row to avoid valgrind warnings in R.Threshold().
2142  R.SetRow(j, infinity());
2143  }
2144  }
2145  R.Threshold(1e-12);
2146 }
2147 
2149  ElementTransformation &Trans,
2150  Vector &dofs) const
2151 {
2152  double data[3];
2153  Vector vk2(data, 2);
2154  Vector vk3(data, 3);
2155 
2156  double * tk_ptr = const_cast<double*>(tk);
2157 
2158  for (int k = 0; k < dof; k++)
2159  {
2160  Trans.SetIntPoint(&Nodes.IntPoint(k));
2161 
2162  vc.Eval(vk3, Trans, Nodes.IntPoint(k));
2163  // dof_k = vk^t J tk
2164  Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
2165  Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2166 
2167  dofs(k) = Trans.Jacobian().InnerProduct(t2, vk2) + t3(2) * vk3(2);
2168  }
2169 
2170 }
2171 
2173  ElementTransformation &Trans,
2174  DenseMatrix &I) const
2175 {
2176  if (fe.GetRangeType() == SCALAR)
2177  {
2178  double vk[Geometry::MaxDim];
2179  Vector shape(fe.GetDof());
2180 
2181  double * tk_ptr = const_cast<double*>(tk);
2182 
2183  I.SetSize(dof, vdim*fe.GetDof());
2184  for (int k = 0; k < dof; k++)
2185  {
2186  const IntegrationPoint &ip = Nodes.IntPoint(k);
2187 
2188  Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
2189  Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2190 
2191  fe.CalcShape(ip, shape);
2192  Trans.SetIntPoint(&ip);
2193  // Transform ND edge tengents from reference to physical space
2194  // vk = J tk
2195  Trans.Jacobian().Mult(t2, vk);
2196  vk[2] = t3[2];
2197  if (fe.GetMapType() == INTEGRAL)
2198  {
2199  double w = 1.0/Trans.Weight();
2200  for (int d = 0; d < vdim; d++)
2201  {
2202  vk[d] *= w;
2203  }
2204  }
2205 
2206  for (int j = 0; j < shape.Size(); j++)
2207  {
2208  double s = shape(j);
2209  if (fabs(s) < 1e-12)
2210  {
2211  s = 0.0;
2212  }
2213  // Project scalar basis function multiplied by each coordinate
2214  // direction onto the transformed edge tangents
2215  for (int d = 0; d < vdim; d++)
2216  {
2217  I(k, j + d*shape.Size()) = s*vk[d];
2218  }
2219  }
2220  }
2221  }
2222  else
2223  {
2224  double vk[Geometry::MaxDim];
2225  DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
2226 
2227  double * tk_ptr = const_cast<double*>(tk);
2228 
2229  I.SetSize(dof, fe.GetDof());
2230  for (int k = 0; k < dof; k++)
2231  {
2232  const IntegrationPoint &ip = Nodes.IntPoint(k);
2233 
2234  Vector t2(&tk_ptr[dof2tk[k] * 3], 2);
2235  Vector t3(&tk_ptr[dof2tk[k] * 3], 3);
2236 
2237  Trans.SetIntPoint(&ip);
2238  // Transform ND edge tangents from reference to physical space
2239  // vk = J tk
2240  Trans.Jacobian().Mult(t2, vk);
2241  // Compute fe basis functions in physical space
2242  fe.CalcVShape(Trans, vshape);
2243  // Project fe basis functions onto transformed edge tangents
2244  for (int j=0; j<vshape.Height(); j++)
2245  {
2246  I(k, j) = 0.0;
2247  for (int i=0; i<2; i++)
2248  {
2249  I(k, j) += vshape(j, i) * vk[i];
2250  }
2251  if (vshape.Width() == 3)
2252  {
2253  I(k, j) += vshape(j, 2) * t3(2);
2254  }
2255  }
2256  }
2257  }
2258 }
2259 
2261  ElementTransformation &Trans,
2262  DenseMatrix &grad) const
2263 {
2264  MFEM_ASSERT(fe.GetMapType() == VALUE, "");
2265 
2266  DenseMatrix dshape(fe.GetDof(), fe.GetDim());
2267  Vector grad_k(fe.GetDof());
2268 
2269  grad.SetSize(dof, fe.GetDof());
2270  for (int k = 0; k < dof; k++)
2271  {
2272  fe.CalcDShape(Nodes.IntPoint(k), dshape);
2273  dshape.Mult(tk + dof2tk[k]*vdim, grad_k);
2274  for (int j = 0; j < grad_k.Size(); j++)
2275  {
2276  grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
2277  }
2278  }
2279 }
2280 
2281 const double ND_R2D_TriangleElement::tk_t[15] =
2282 { 1.,0.,0., -1.,1.,0., 0.,-1.,0., 0.,1.,0., 0.,0.,1. };
2283 
2285  const int cb_type)
2286  : ND_R2D_FiniteElement(p, Geometry::TRIANGLE, ((3*p + 1)*(p + 2))/2, tk_t),
2287  ND_FE(p),
2288  H1_FE(p, cb_type)
2289 {
2290  int pm1 = p - 1, pm2 = p - 2;
2291 
2292 #ifndef MFEM_THREAD_SAFE
2293  nd_shape.SetSize(ND_FE.GetDof(), 2);
2294  h1_shape.SetSize(H1_FE.GetDof());
2295  nd_dshape.SetSize(ND_FE.GetDof(), 1);
2296  h1_dshape.SetSize(H1_FE.GetDof(), 2);
2297 #endif
2298 
2299  int o = 0;
2300  int n = 0;
2301  int h = 0;
2302  // Three nodes
2303  dof_map[o] = -1 - h++; dof2tk[o++] = 4;
2304  dof_map[o] = -1 - h++; dof2tk[o++] = 4;
2305  dof_map[o] = -1 - h++; dof2tk[o++] = 4;
2306 
2307  // Three edges
2308  for (int e=0; e<3; e++)
2309  {
2310  // Dofs in the plane
2311  for (int i=0; i<p; i++)
2312  {
2313  dof_map[o] = n++; dof2tk[o++] = e;
2314  }
2315  // z-directed dofs
2316  for (int i=0; i<pm1; i++)
2317  {
2318  dof_map[o] = -1 - h++; dof2tk[o++] = 4;
2319  }
2320  }
2321 
2322  // Interior dofs in the plane
2323  for (int j = 0; j <= pm2; j++)
2324  for (int i = 0; i + j <= pm2; i++)
2325  {
2326  dof_map[o] = n++; dof2tk[o++] = 0;
2327  dof_map[o] = n++; dof2tk[o++] = 3;
2328  }
2329 
2330  // Interior z-directed dofs
2331  for (int j = 0; j < pm1; j++)
2332  for (int i = 0; i + j < pm2; i++)
2333  {
2334  dof_map[o] = -1 - h++; dof2tk[o++] = 4;
2335  }
2336 
2337  MFEM_VERIFY(n == ND_FE.GetDof(),
2338  "ND_R2D_Triangle incorrect number of ND dofs.");
2339  MFEM_VERIFY(h == H1_FE.GetDof(),
2340  "ND_R2D_Triangle incorrect number of H1 dofs.");
2341  MFEM_VERIFY(o == GetDof(),
2342  "ND_R2D_Triangle incorrect number of dofs.");
2343 
2344  const IntegrationRule & nd_Nodes = ND_FE.GetNodes();
2345  const IntegrationRule & h1_Nodes = H1_FE.GetNodes();
2346 
2347  for (int i=0; i<dof; i++)
2348  {
2349  int idx = dof_map[i];
2350  if (idx >= 0)
2351  {
2352  const IntegrationPoint & ip = nd_Nodes.IntPoint(idx);
2353  Nodes.IntPoint(i).Set2(ip.x, ip.y);
2354  }
2355  else
2356  {
2357  const IntegrationPoint & ip = h1_Nodes.IntPoint(-idx-1);
2358  Nodes.IntPoint(i).Set2(ip.x, ip.y);
2359  }
2360  }
2361 }
2362 
2364  DenseMatrix &shape) const
2365 {
2366 #ifdef MFEM_THREAD_SAFE
2367  DenseMatrix nd_shape(ND_FE.GetDof(), 2);
2368  Vector h1_shape(H1_FE.GetDof());
2369 #endif
2370 
2371  ND_FE.CalcVShape(ip, nd_shape);
2372  H1_FE.CalcShape(ip, h1_shape);
2373 
2374  for (int i=0; i<dof; i++)
2375  {
2376  int idx = dof_map[i];
2377  if (idx >= 0)
2378  {
2379  shape(i, 0) = nd_shape(idx, 0);
2380  shape(i, 1) = nd_shape(idx, 1);
2381  shape(i, 2) = 0.0;
2382  }
2383  else
2384  {
2385  shape(i, 0) = 0.0;
2386  shape(i, 1) = 0.0;
2387  shape(i, 2) = h1_shape(-idx-1);
2388  }
2389  }
2390 }
2391 
2393  DenseMatrix &curl_shape) const
2394 {
2395 #ifdef MFEM_THREAD_SAFE
2396  DenseMatrix nd_dshape(ND_FE.GetDof(), 1);
2397  DenseMatrix h1_dshape(H1_FE.GetDof(), 2);
2398 #endif
2399 
2400  ND_FE.CalcCurlShape(ip, nd_dshape);
2401  H1_FE.CalcDShape(ip, h1_dshape);
2402 
2403  for (int i=0; i<dof; i++)
2404  {
2405  int idx = dof_map[i];
2406  if (idx >= 0)
2407  {
2408  curl_shape(i, 0) = 0.0;
2409  curl_shape(i, 1) = 0.0;
2410  curl_shape(i, 2) = nd_dshape(idx, 0);
2411  }
2412  else
2413  {
2414  curl_shape(i, 0) = h1_dshape(-idx-1, 1);
2415  curl_shape(i, 1) = -h1_dshape(-idx-1, 0);
2416  curl_shape(i, 2) = 0.0;
2417  }
2418  }
2419 }
2420 
2421 
2422 const double ND_R2D_QuadrilateralElement::tk_q[15] =
2423 { 1.,0.,0., 0.,1.,0., -1.,0.,0., 0.,-1.,0., 0.,0.,1. };
2424 
2426  const int cb_type,
2427  const int ob_type)
2428  : ND_R2D_FiniteElement(p, Geometry::SQUARE, ((3*p + 1)*(p + 1)), tk_q),
2429  cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
2430  obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
2431 {
2432  const double *cp = poly1d.ClosedPoints(p, cb_type);
2433  const double *op = poly1d.OpenPoints(p - 1, ob_type);
2434  const int dofx = p*(p+1);
2435  const int dofy = p*(p+1);
2436  const int dofxy = dofx+dofy;
2437 
2438 #ifndef MFEM_THREAD_SAFE
2439  shape_cx.SetSize(p + 1);
2440  shape_ox.SetSize(p);
2441  shape_cy.SetSize(p + 1);
2442  shape_oy.SetSize(p);
2443  dshape_cx.SetSize(p + 1);
2444  dshape_cy.SetSize(p + 1);
2445 #endif
2446 
2447  dof_map.SetSize(dof);
2448 
2449  int o = 0;
2450  // nodes
2451  dof_map[dofxy] = o++; // (0)
2452  dof_map[dofxy+p] = o++; // (1)
2453  dof_map[dof-1] = o++; // (2)
2454  dof_map[dof-p-1] = o++; // (3)
2455 
2456  // edges
2457  for (int i = 0; i < p; i++) // (0,1) - x-directed
2458  {
2459  dof_map[i + 0*p] = o++;
2460  }
2461  for (int i = 1; i < p; i++) // (0,1) - z-directed
2462  {
2463  dof_map[dofxy + i + 0*(p+1)] = o++;
2464  }
2465  for (int j = 0; j < p; j++) // (1,2) - y-directed
2466  {
2467  dof_map[dofx + p + j*(p + 1)] = o++;
2468  }
2469  for (int j = 1; j < p; j++) // (1,2) - z-directed
2470  {
2471  dof_map[dofxy + p + j*(p + 1)] = o++;
2472  }
2473  for (int i = 0; i < p; i++) // (2,3) - x-directed
2474  {
2475  dof_map[(p - 1 - i) + p*p] = -1 - (o++);
2476  }
2477  for (int i = 1; i < p; i++) // (2,3) - z-directed
2478  {
2479  dof_map[dofxy + (p - i) + p*(p + 1)] = o++;
2480  }
2481  for (int j = 0; j < p; j++) // (3,0) - y-directed
2482  {
2483  dof_map[dofx + 0 + (p - 1 - j)*(p + 1)] = -1 - (o++);
2484  }
2485  for (int j = 1; j < p; j++) // (3,0) - z-directed
2486  {
2487  dof_map[dofxy + (p - j)*(p + 1)] = o++;
2488  }
2489 
2490  // interior
2491  // x-components
2492  for (int j = 1; j < p; j++)
2493  for (int i = 0; i < p; i++)
2494  {
2495  dof_map[i + j*p] = o++;
2496  }
2497  // y-components
2498  for (int j = 0; j < p; j++)
2499  for (int i = 1; i < p; i++)
2500  {
2501  dof_map[dofx + i + j*(p + 1)] = o++;
2502  }
2503  // z-components
2504  for (int j = 1; j < p; j++)
2505  for (int i = 1; i < p; i++)
2506  {
2507  dof_map[dofxy + i + j*(p + 1)] = o++;
2508  }
2509 
2510  // set dof2tk and Nodes
2511  o = 0;
2512  // x-components
2513  for (int j = 0; j <= p; j++)
2514  for (int i = 0; i < p; i++)
2515  {
2516  int idx;
2517  if ((idx = dof_map[o++]) < 0)
2518  {
2519  dof2tk[idx = -1 - idx] = 2;
2520  }
2521  else
2522  {
2523  dof2tk[idx] = 0;
2524  }
2525  Nodes.IntPoint(idx).Set2(op[i], cp[j]);
2526  }
2527  // y-components
2528  for (int j = 0; j < p; j++)
2529  for (int i = 0; i <= p; i++)
2530  {
2531  int idx;
2532  if ((idx = dof_map[o++]) < 0)
2533  {
2534  dof2tk[idx = -1 - idx] = 3;
2535  }
2536  else
2537  {
2538  dof2tk[idx] = 1;
2539  }
2540  Nodes.IntPoint(idx).Set2(cp[i], op[j]);
2541  }
2542  // z-components
2543  for (int j = 0; j <= p; j++)
2544  for (int i = 0; i <= p; i++)
2545  {
2546  int idx = dof_map[o++];
2547  dof2tk[idx] = 4;
2548  Nodes.IntPoint(idx).Set2(cp[i], cp[j]);
2549  }
2550 }
2551 
2553  DenseMatrix &shape) const
2554 {
2555  const int p = order;
2556 
2557 #ifdef MFEM_THREAD_SAFE
2558  Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
2559 #endif
2560 
2561  cbasis1d.Eval(ip.x, shape_cx);
2562  obasis1d.Eval(ip.x, shape_ox);
2563  cbasis1d.Eval(ip.y, shape_cy);
2564  obasis1d.Eval(ip.y, shape_oy);
2565 
2566  int o = 0;
2567  // x-components
2568  for (int j = 0; j <= p; j++)
2569  for (int i = 0; i < p; i++)
2570  {
2571  int idx, s;
2572  if ((idx = dof_map[o++]) < 0)
2573  {
2574  idx = -1 - idx, s = -1;
2575  }
2576  else
2577  {
2578  s = +1;
2579  }
2580  shape(idx,0) = s*shape_ox(i)*shape_cy(j);
2581  shape(idx,1) = 0.;
2582  shape(idx,2) = 0.;
2583  }
2584  // y-components
2585  for (int j = 0; j < p; j++)
2586  for (int i = 0; i <= p; i++)
2587  {
2588  int idx, s;
2589  if ((idx = dof_map[o++]) < 0)
2590  {
2591  idx = -1 - idx, s = -1;
2592  }
2593  else
2594  {
2595  s = +1;
2596  }
2597  shape(idx,0) = 0.;
2598  shape(idx,1) = s*shape_cx(i)*shape_oy(j);
2599  shape(idx,2) = 0.;
2600  }
2601  // z-components
2602  for (int j = 0; j <= p; j++)
2603  for (int i = 0; i <= p; i++)
2604  {
2605  int idx = dof_map[o++];
2606  shape(idx,0) = 0.;
2607  shape(idx,1) = 0.;
2608  shape(idx,2) = shape_cx(i)*shape_cy(j);
2609  }
2610 }
2611 
2613  DenseMatrix &curl_shape) const
2614 {
2615  const int p = order;
2616 
2617 #ifdef MFEM_THREAD_SAFE
2618  Vector shape_cx(p + 1), shape_ox(p), shape_cy(p + 1), shape_oy(p);
2619  Vector dshape_cx(p + 1), dshape_cy(p + 1);
2620 #endif
2621 
2622  cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
2623  obasis1d.Eval(ip.x, shape_ox);
2624  cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
2625  obasis1d.Eval(ip.y, shape_oy);
2626 
2627  int o = 0;
2628  // x-components
2629  for (int j = 0; j <= p; j++)
2630  for (int i = 0; i < p; i++)
2631  {
2632  int idx, s;
2633  if ((idx = dof_map[o++]) < 0)
2634  {
2635  idx = -1 - idx, s = -1;
2636  }
2637  else
2638  {
2639  s = +1;
2640  }
2641  curl_shape(idx,0) = 0.;
2642  curl_shape(idx,1) = 0.;
2643  curl_shape(idx,2) = -s*shape_ox(i)*dshape_cy(j);
2644  }
2645  // y-components
2646  for (int j = 0; j < p; j++)
2647  for (int i = 0; i <= p; i++)
2648  {
2649  int idx, s;
2650  if ((idx = dof_map[o++]) < 0)
2651  {
2652  idx = -1 - idx, s = -1;
2653  }
2654  else
2655  {
2656  s = +1;
2657  }
2658  curl_shape(idx,0) = 0.;
2659  curl_shape(idx,1) = 0.;
2660  curl_shape(idx,2) = s*dshape_cx(i)*shape_oy(j);
2661  }
2662  // z-components
2663  for (int j = 0; j <= p; j++)
2664  for (int i = 0; i <= p; i++)
2665  {
2666  int idx = dof_map[o++];
2667  curl_shape(idx,0) = shape_cx(i)*dshape_cy(j);
2668  curl_shape(idx,1) = -dshape_cx(i)*shape_cy(j);
2669  curl_shape(idx,2) = 0.;
2670  }
2671 }
2672 
2673 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
void Set(const double *p, const int dim)
Definition: intrules.hpp:37
ND_SegmentElement(const int p, const int ob_type=BasisType::GaussLegendre)
Construct the ND_SegmentElement of order p and open BasisType ob_type.
Definition: fe_nd.cpp:1198
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition: fe_nd.cpp:2039
void ProjectIntegrated(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Definition: fe_nd.cpp:576
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_base.cpp:39
void ScaleIntegrated(bool scale_integrated_)
Set whether the &quot;integrated&quot; basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition: fe_base.cpp:1879
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
bool IsIntegratedType() const
Returns true if the basis is &quot;integrated&quot;, false otherwise.
Definition: fe_base.hpp:1022
Base class for vector Coefficients that optionally depend on time and space.
void SetRow(int r, const double *row)
Definition: densemat.cpp:2139
ND_R1D_PointElement(int p)
Construct the ND_R1D_PointElement.
Definition: fe_nd.cpp:1513
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
int dim
Dimension of reference space.
Definition: fe_base.hpp:238
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:297
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
ND_R2D_TriangleElement(const int p, const int cb_type=BasisType::GaussLobatto)
Construct the ND_R2D_TriangleElement of order p.
Definition: fe_nd.cpp:2284
ND_R2D_SegmentElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R2D_SegmentElement of order p and closed and open BasisType cb_type and ob_type...
Definition: fe_nd.cpp:1819
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:1868
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:145
virtual void CalcVShape(ElementTransformation &Trans, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in physical space at the point ...
Definition: fe_nd.cpp:2022
DenseMatrix vshape
Definition: fe_base.hpp:250
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:349
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3878
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition: fe_base.hpp:241
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition: eltrans.hpp:98
Intermediate class for finite elements whose basis functions return vector values.
Definition: fe_base.hpp:786
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:71
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:2363
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:1606
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:1911
const double * ClosedPoints(const int p, const int btype=BasisType::GaussLobatto)
Get coordinates of a closed (GaussLegendre) set of points if degree p.
Definition: fe_base.hpp:1070
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:1524
ND_WedgeElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Definition: fe_nd.cpp:1226
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:1659
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:1119
ND_TriangleElement(const int p)
Construct the ND_TriangleElement of order p.
Definition: fe_nd.cpp:1038
int cdim
Dimension of curl for vector-valued basis functions.
Definition: fe_base.hpp:240
int vdim
Vector dimension of vector-valued basis functions.
Definition: fe_base.hpp:239
static void CalcBasis(const int p, const double x, double *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition: fe_base.hpp:1086
ND_R2D_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R2D_QuadrilateralElement of order p and closed and open BasisType cb_type and ob_typ...
Definition: fe_nd.cpp:2425
Geometry Geometries
Definition: fe.cpp:49
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
ND_HexahedronElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_HexahedronElement of order p and closed and open BasisType cb_type and ob_type...
Definition: fe_nd.cpp:25
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:533
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition: fe_nd.cpp:1699
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
Definition: fe_nd.cpp:2148
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
ND_R2D_FiniteElement(int p, Geometry::Type G, int Do, const double *tk_fe)
Definition: fe_nd.cpp:2004
ND_TetrahedronElement(const int p)
Construct the ND_TetrahedronElement of order p.
Definition: fe_nd.cpp:776
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
const double * tk
Definition: fe_nd.hpp:552
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the &quot;integrated&quot; basis type using pre-computed closed basis derivatives. ...
Definition: fe_base.cpp:1854
ND_R1D_SegmentElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_R1D_SegmentElement of order p and closed and open BasisType cb_type and ob_type...
Definition: fe_nd.cpp:1542
void Set2(const double x1, const double x2)
Definition: intrules.hpp:80
int GetVDim()
Returns dimension of the vector.
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:2392
const double * OpenPoints(const int p, const int btype=BasisType::GaussLegendre)
Get coordinates of an open (GaussLegendre) set of points if degree p.
Definition: fe_base.hpp:1065
static const int MaxDim
Definition: geom.hpp:43
void Threshold(double eps)
Replace small entries, abs(a_ij) &lt;= eps, with zero.
Definition: densemat.cpp:2169
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
Definition: fe_base.hpp:314
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:713
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
IntegrationRule Nodes
Definition: fe_base.hpp:248
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:1153
ND_QuadrilateralElement(const int p, const int cb_type=BasisType::GaussLobatto, const int ob_type=BasisType::GaussLegendre)
Construct the ND_QuadrilateralElement of order p and closed and open BasisType cb_type and ob_type...
Definition: fe_nd.cpp:485
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
Definition: fe_nd.cpp:1710
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3924
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:309
virtual void Project(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector coefficient and a transformation, compute its projection (approximation) in the local ...
Definition: fe_nd.cpp:1979
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:1475
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:2612
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_base.cpp:64
Class for integration point with weight.
Definition: intrules.hpp:25
int dof
Number of degrees of freedom.
Definition: fe_base.hpp:245
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:924
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:652
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:2552
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:1441
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:59
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
RefCoord t[3]
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs. ...
Definition: fe_nd.cpp:2102
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
Definition: fe_h1.cpp:40
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_nd.cpp:1215
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:397
void Eval(const double x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition: fe_base.cpp:1622
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
Definition: fe_base.cpp:611
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
virtual void ProjectGrad(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_nd.cpp:2260
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:405
Describes the function space on each element.
Definition: fe_base.hpp:218
RefCoord s[3]
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_nd.cpp:966
Array< int > dof_map
Definition: fe_nd.hpp:553
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe_base.hpp:340
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
void ProjectIntegrated(VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Definition: fe_nd.cpp:244
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:245
Poly_1D poly1d
Definition: fe.cpp:28
Implements CalcCurlShape methods.
Definition: fe_base.hpp:297
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_h1.cpp:555