MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fe_pos.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 // H1 Finite Element classes utilizing the Bernstein basis
13 
14 #include "fe_pos.hpp"
15 #include "../bilininteg.hpp"
16 #include "../coefficient.hpp"
17 
18 namespace mfem
19 {
20 
21 using namespace std;
22 
24  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
25 {
26  for (int i = 0; i < dof; i++)
27  {
28  const IntegrationPoint &ip = Nodes.IntPoint(i);
29  Trans.SetIntPoint(&ip);
30  dofs(i) = coeff.Eval(Trans, ip);
31  }
32 }
33 
36 {
37  MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
38  Vector x(vc.GetVDim());
39 
40  for (int i = 0; i < dof; i++)
41  {
42  const IntegrationPoint &ip = Nodes.IntPoint(i);
43  Trans.SetIntPoint(&ip);
44  vc.Eval (x, Trans, ip);
45  for (int j = 0; j < x.Size(); j++)
46  {
47  dofs(dof*j+i) = x(j);
48  }
49  }
50 }
51 
54 {
55  const NodalFiniteElement *nfe =
56  dynamic_cast<const NodalFiniteElement *>(&fe);
57 
58  if (nfe && dof == nfe->GetDof())
59  {
60  nfe->Project(*this, Trans, I);
61  I.Invert();
62  }
63  else
64  {
65  // local L2 projection
66  DenseMatrix pos_mass, mixed_mass;
67  MassIntegrator mass_integ;
68 
69  mass_integ.AssembleElementMatrix(*this, Trans, pos_mass);
70  mass_integ.AssembleElementMatrix2(fe, *this, Trans, mixed_mass);
71 
72  DenseMatrixInverse pos_mass_inv(pos_mass);
73  I.SetSize(dof, fe.GetDof());
74  pos_mass_inv.Mult(mixed_mass, I);
75  }
76 }
77 
78 
80  const int dims, const int p, const DofMapType dmtype)
81  : PositiveFiniteElement(dims, GetTensorProductGeometry(dims),
82  Pow(p + 1, dims), p,
83  dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk),
84  TensorBasisElement(dims, p, BasisType::Positive, dmtype) { }
85 
86 
88  : PositiveFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
89 {
90  Nodes.IntPoint(0).x = 0.0;
91  Nodes.IntPoint(0).y = 0.0;
92  Nodes.IntPoint(1).x = 1.0;
93  Nodes.IntPoint(1).y = 0.0;
94  Nodes.IntPoint(2).x = 1.0;
95  Nodes.IntPoint(2).y = 1.0;
96  Nodes.IntPoint(3).x = 0.0;
97  Nodes.IntPoint(3).y = 1.0;
98  Nodes.IntPoint(4).x = 0.5;
99  Nodes.IntPoint(4).y = 0.0;
100  Nodes.IntPoint(5).x = 1.0;
101  Nodes.IntPoint(5).y = 0.5;
102  Nodes.IntPoint(6).x = 0.5;
103  Nodes.IntPoint(6).y = 1.0;
104  Nodes.IntPoint(7).x = 0.0;
105  Nodes.IntPoint(7).y = 0.5;
106  Nodes.IntPoint(8).x = 0.5;
107  Nodes.IntPoint(8).y = 0.5;
108 }
109 
111  Vector &shape) const
112 {
113  double x = ip.x, y = ip.y;
114  double l1x, l2x, l3x, l1y, l2y, l3y;
115 
116  l1x = (1. - x) * (1. - x);
117  l2x = 2. * x * (1. - x);
118  l3x = x * x;
119  l1y = (1. - y) * (1. - y);
120  l2y = 2. * y * (1. - y);
121  l3y = y * y;
122 
123  shape(0) = l1x * l1y;
124  shape(4) = l2x * l1y;
125  shape(1) = l3x * l1y;
126  shape(7) = l1x * l2y;
127  shape(8) = l2x * l2y;
128  shape(5) = l3x * l2y;
129  shape(3) = l1x * l3y;
130  shape(6) = l2x * l3y;
131  shape(2) = l3x * l3y;
132 }
133 
135  DenseMatrix &dshape) const
136 {
137  double x = ip.x, y = ip.y;
138  double l1x, l2x, l3x, l1y, l2y, l3y;
139  double d1x, d2x, d3x, d1y, d2y, d3y;
140 
141  l1x = (1. - x) * (1. - x);
142  l2x = 2. * x * (1. - x);
143  l3x = x * x;
144  l1y = (1. - y) * (1. - y);
145  l2y = 2. * y * (1. - y);
146  l3y = y * y;
147 
148  d1x = 2. * x - 2.;
149  d2x = 2. - 4. * x;
150  d3x = 2. * x;
151  d1y = 2. * y - 2.;
152  d2y = 2. - 4. * y;
153  d3y = 2. * y;
154 
155  dshape(0,0) = d1x * l1y;
156  dshape(0,1) = l1x * d1y;
157 
158  dshape(4,0) = d2x * l1y;
159  dshape(4,1) = l2x * d1y;
160 
161  dshape(1,0) = d3x * l1y;
162  dshape(1,1) = l3x * d1y;
163 
164  dshape(7,0) = d1x * l2y;
165  dshape(7,1) = l1x * d2y;
166 
167  dshape(8,0) = d2x * l2y;
168  dshape(8,1) = l2x * d2y;
169 
170  dshape(5,0) = d3x * l2y;
171  dshape(5,1) = l3x * d2y;
172 
173  dshape(3,0) = d1x * l3y;
174  dshape(3,1) = l1x * d3y;
175 
176  dshape(6,0) = d2x * l3y;
177  dshape(6,1) = l2x * d3y;
178 
179  dshape(2,0) = d3x * l3y;
180  dshape(2,1) = l3x * d3y;
181 }
182 
185 {
186  double s[9];
187  IntegrationPoint tr_ip;
188  Vector xx(&tr_ip.x, 2), shape(s, 9);
189 
190  for (int i = 0; i < 9; i++)
191  {
192  Trans.Transform(Nodes.IntPoint(i), xx);
193  CalcShape(tr_ip, shape);
194  for (int j = 0; j < 9; j++)
195  if (fabs(I(i,j) = s[j]) < 1.0e-12)
196  {
197  I(i,j) = 0.0;
198  }
199  }
200  for (int i = 0; i < 9; i++)
201  {
202  double *d = &I(0,i);
203  d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
204  d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
205  d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
206  d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
207  d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
208  0.25 * (d[0] + d[1] + d[2] + d[3]);
209  }
210 }
211 
213  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
214 {
215  double *d = dofs;
216 
217  for (int i = 0; i < 9; i++)
218  {
219  const IntegrationPoint &ip = Nodes.IntPoint(i);
220  Trans.SetIntPoint(&ip);
221  d[i] = coeff.Eval(Trans, ip);
222  }
223  d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
224  d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
225  d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
226  d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
227  d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
228  0.25 * (d[0] + d[1] + d[2] + d[3]);
229 }
230 
233  Vector &dofs) const
234 {
235  double v[3];
236  Vector x (v, vc.GetVDim());
237 
238  for (int i = 0; i < 9; i++)
239  {
240  const IntegrationPoint &ip = Nodes.IntPoint(i);
241  Trans.SetIntPoint(&ip);
242  vc.Eval (x, Trans, ip);
243  for (int j = 0; j < x.Size(); j++)
244  {
245  dofs(9*j+i) = v[j];
246  }
247  }
248  for (int j = 0; j < x.Size(); j++)
249  {
250  double *d = &dofs(9*j);
251 
252  d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
253  d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
254  d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
255  d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
256  d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
257  0.25 * (d[0] + d[1] + d[2] + d[3]);
258  }
259 }
260 
261 
263  : PositiveFiniteElement(1, Geometry::SEGMENT, 3, 2)
264 {
265  Nodes.IntPoint(0).x = 0.0;
266  Nodes.IntPoint(1).x = 1.0;
267  Nodes.IntPoint(2).x = 0.5;
268 }
269 
271  Vector &shape) const
272 {
273  const double x = ip.x, x1 = 1. - x;
274 
275  shape(0) = x1 * x1;
276  shape(1) = x * x;
277  shape(2) = 2. * x * x1;
278 }
279 
281  DenseMatrix &dshape) const
282 {
283  const double x = ip.x;
284 
285  dshape(0,0) = 2. * x - 2.;
286  dshape(1,0) = 2. * x;
287  dshape(2,0) = 2. - 4. * x;
288 }
289 
290 
292  : PositiveTensorFiniteElement(1, p, H1_DOF_MAP)
293 {
294 #ifndef MFEM_THREAD_SAFE
295  // thread private versions; see class header.
296  shape_x.SetSize(p+1);
297  dshape_x.SetSize(p+1);
298 #endif
299 
300  // Endpoints need to be first in the list, so reorder them.
301  Nodes.IntPoint(0).x = 0.0;
302  Nodes.IntPoint(1).x = 1.0;
303  for (int i = 1; i < p; i++)
304  {
305  Nodes.IntPoint(i+1).x = double(i)/p;
306  }
307 }
308 
310  Vector &shape) const
311 {
312  const int p = order;
313 
314 #ifdef MFEM_THREAD_SAFE
315  Vector shape_x(p+1);
316 #endif
317 
318  Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData() );
319 
320  // Endpoints need to be first in the list, so reorder them.
321  shape(0) = shape_x(0);
322  shape(1) = shape_x(p);
323  for (int i = 1; i < p; i++)
324  {
325  shape(i+1) = shape_x(i);
326  }
327 }
328 
330  DenseMatrix &dshape) const
331 {
332  const int p = order;
333 
334 #ifdef MFEM_THREAD_SAFE
335  Vector shape_x(p+1), dshape_x(p+1);
336 #endif
337 
338  Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData(), dshape_x.GetData() );
339 
340  // Endpoints need to be first in the list, so reorder them.
341  dshape(0,0) = dshape_x(0);
342  dshape(1,0) = dshape_x(p);
343  for (int i = 1; i < p; i++)
344  {
345  dshape(i+1,0) = dshape_x(i);
346  }
347 }
348 
349 void H1Pos_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
350 {
351  dofs = 0.0;
352  dofs[vertex] = 1.0;
353 }
354 
355 
357  : PositiveTensorFiniteElement(2, p, H1_DOF_MAP)
358 {
359 #ifndef MFEM_THREAD_SAFE
360  const int p1 = p + 1;
361 
362  shape_x.SetSize(p1);
363  shape_y.SetSize(p1);
364  dshape_x.SetSize(p1);
365  dshape_y.SetSize(p1);
366 #endif
367 
368  int o = 0;
369  for (int j = 0; j <= p; j++)
370  for (int i = 0; i <= p; i++)
371  {
372  Nodes.IntPoint(dof_map[o++]).Set2(double(i)/p, double(j)/p);
373  }
374 }
375 
377  Vector &shape) const
378 {
379  const int p = order;
380 
381 #ifdef MFEM_THREAD_SAFE
382  Vector shape_x(p+1), shape_y(p+1);
383 #endif
384 
385  Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData() );
386  Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData() );
387 
388  // Reorder so that vertices are at the beginning of the list
389  for (int o = 0, j = 0; j <= p; j++)
390  for (int i = 0; i <= p; i++)
391  {
392  shape(dof_map[o++]) = shape_x(i)*shape_y(j);
393  }
394 }
395 
397  DenseMatrix &dshape) const
398 {
399  const int p = order;
400 
401 #ifdef MFEM_THREAD_SAFE
402  Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
403 #endif
404 
405  Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData(), dshape_x.GetData() );
406  Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData(), dshape_y.GetData() );
407 
408  // Reorder so that vertices are at the beginning of the list
409  for (int o = 0, j = 0; j <= p; j++)
410  for (int i = 0; i <= p; i++)
411  {
412  dshape(dof_map[o],0) = dshape_x(i)* shape_y(j);
413  dshape(dof_map[o],1) = shape_x(i)*dshape_y(j); o++;
414  }
415 }
416 
418 {
419  dofs = 0.0;
420  dofs[vertex] = 1.0;
421 }
422 
423 
425  : PositiveTensorFiniteElement(3, p, H1_DOF_MAP)
426 {
427 #ifndef MFEM_THREAD_SAFE
428  const int p1 = p + 1;
429 
430  shape_x.SetSize(p1);
431  shape_y.SetSize(p1);
432  shape_z.SetSize(p1);
433  dshape_x.SetSize(p1);
434  dshape_y.SetSize(p1);
435  dshape_z.SetSize(p1);
436 #endif
437 
438  int o = 0;
439  for (int k = 0; k <= p; k++)
440  for (int j = 0; j <= p; j++)
441  for (int i = 0; i <= p; i++)
442  Nodes.IntPoint(dof_map[o++]).Set3(double(i)/p, double(j)/p,
443  double(k)/p);
444 }
445 
447  Vector &shape) const
448 {
449  const int p = order;
450 
451 #ifdef MFEM_THREAD_SAFE
452  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
453 #endif
454 
455  Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData() );
456  Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData() );
457  Poly_1D::CalcBernstein(p, ip.z, shape_z.GetData() );
458 
459  for (int o = 0, k = 0; k <= p; k++)
460  for (int j = 0; j <= p; j++)
461  for (int i = 0; i <= p; i++)
462  {
463  shape(dof_map[o++]) = shape_x(i)*shape_y(j)*shape_z(k);
464  }
465 }
466 
468  DenseMatrix &dshape) const
469 {
470  const int p = order;
471 
472 #ifdef MFEM_THREAD_SAFE
473  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
474  Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
475 #endif
476 
477  Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData(), dshape_x.GetData() );
478  Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData(), dshape_y.GetData() );
479  Poly_1D::CalcBernstein(p, ip.z, shape_z.GetData(), dshape_z.GetData() );
480 
481  for (int o = 0, k = 0; k <= p; k++)
482  for (int j = 0; j <= p; j++)
483  for (int i = 0; i <= p; i++)
484  {
485  dshape(dof_map[o],0) = dshape_x(i)* shape_y(j)* shape_z(k);
486  dshape(dof_map[o],1) = shape_x(i)*dshape_y(j)* shape_z(k);
487  dshape(dof_map[o],2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
488  }
489 }
490 
491 void H1Pos_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
492 {
493  dofs = 0.0;
494  dofs[vertex] = 1.0;
495 }
496 
497 
499  : PositiveFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
500  FunctionSpace::Pk)
501 {
502 #ifndef MFEM_THREAD_SAFE
504  dshape_1d.SetSize(p + 1);
506 #endif
508 
509  struct Index
510  {
511  int p2p3;
512  Index(int p) { p2p3 = 2*p + 3; }
513  int operator()(int i, int j) { return ((p2p3-j)*j)/2+i; }
514  };
515  Index idx(p);
516 
517  // vertices
518  dof_map[idx(0,0)] = 0;
519  Nodes.IntPoint(0).Set2(0., 0.);
520  dof_map[idx(p,0)] = 1;
521  Nodes.IntPoint(1).Set2(1., 0.);
522  dof_map[idx(0,p)] = 2;
523  Nodes.IntPoint(2).Set2(0., 1.);
524 
525  // edges
526  int o = 3;
527  for (int i = 1; i < p; i++)
528  {
529  dof_map[idx(i,0)] = o;
530  Nodes.IntPoint(o++).Set2(double(i)/p, 0.);
531  }
532  for (int i = 1; i < p; i++)
533  {
534  dof_map[idx(p-i,i)] = o;
535  Nodes.IntPoint(o++).Set2(double(p-i)/p, double(i)/p);
536  }
537  for (int i = 1; i < p; i++)
538  {
539  dof_map[idx(0,p-i)] = o;
540  Nodes.IntPoint(o++).Set2(0., double(p-i)/p);
541  }
542 
543  // interior
544  for (int j = 1; j < p; j++)
545  for (int i = 1; i + j < p; i++)
546  {
547  dof_map[idx(i,j)] = o;
548  Nodes.IntPoint(o++).Set2(double(i)/p, double(j)/p);
549  }
550 }
551 
552 // static method
554  const int p, const double l1, const double l2, double *shape)
555 {
556  const double l3 = 1. - l1 - l2;
557 
558  // The (i,j) basis function is given by: T(i,j,p-i-j) l1^i l2^j l3^{p-i-j},
559  // where T(i,j,k) = (i+j+k)! / (i! j! k!)
560  // Another expression is given by the terms of the expansion:
561  // (l1 + l2 + l3)^p =
562  // \sum_{j=0}^p \binom{p}{j} l2^j
563  // \sum_{i=0}^{p-j} \binom{p-j}{i} l1^i l3^{p-j-i}
564  const int *bp = Poly_1D::Binom(p);
565  double z = 1.;
566  for (int o = 0, j = 0; j <= p; j++)
567  {
568  Poly_1D::CalcBinomTerms(p - j, l1, l3, &shape[o]);
569  double s = bp[j]*z;
570  for (int i = 0; i <= p - j; i++)
571  {
572  shape[o++] *= s;
573  }
574  z *= l2;
575  }
576 }
577 
578 // static method
580  const int p, const double l1, const double l2,
581  double *dshape_1d, double *dshape)
582 {
583  const int dof = ((p + 1)*(p + 2))/2;
584  const double l3 = 1. - l1 - l2;
585 
586  const int *bp = Poly_1D::Binom(p);
587  double z = 1.;
588  for (int o = 0, j = 0; j <= p; j++)
589  {
590  Poly_1D::CalcDBinomTerms(p - j, l1, l3, dshape_1d);
591  double s = bp[j]*z;
592  for (int i = 0; i <= p - j; i++)
593  {
594  dshape[o++] = s*dshape_1d[i];
595  }
596  z *= l2;
597  }
598  z = 1.;
599  for (int i = 0; i <= p; i++)
600  {
601  Poly_1D::CalcDBinomTerms(p - i, l2, l3, dshape_1d);
602  double s = bp[i]*z;
603  for (int o = i, j = 0; j <= p - i; j++)
604  {
605  dshape[dof + o] = s*dshape_1d[j];
606  o += p + 1 - j;
607  }
608  z *= l1;
609  }
610 }
611 
613  Vector &shape) const
614 {
615 #ifdef MFEM_THREAD_SAFE
616  Vector m_shape(dof);
617 #endif
618  CalcShape(order, ip.x, ip.y, m_shape.GetData());
619  for (int i = 0; i < dof; i++)
620  {
621  shape(dof_map[i]) = m_shape(i);
622  }
623 }
624 
626  DenseMatrix &dshape) const
627 {
628 #ifdef MFEM_THREAD_SAFE
629  Vector dshape_1d(order + 1);
631 #endif
632  CalcDShape(order, ip.x, ip.y, dshape_1d.GetData(), m_dshape.Data());
633  for (int d = 0; d < 2; d++)
634  {
635  for (int i = 0; i < dof; i++)
636  {
637  dshape(dof_map[i],d) = m_dshape(i,d);
638  }
639  }
640 }
641 
642 
644  : PositiveFiniteElement(3, Geometry::TETRAHEDRON,
645  ((p + 1)*(p + 2)*(p + 3))/6, p, FunctionSpace::Pk)
646 {
647 #ifndef MFEM_THREAD_SAFE
649  dshape_1d.SetSize(p + 1);
651 #endif
653 
654  struct Index
655  {
656  int p, dof;
657  int tri(int k) { return (k*(k + 1))/2; }
658  int tet(int k) { return (k*(k + 1)*(k + 2))/6; }
659  Index(int p_) { p = p_; dof = tet(p + 1); }
660  int operator()(int i, int j, int k)
661  { return dof - tet(p - k) - tri(p + 1 - k - j) + i; }
662  };
663  Index idx(p);
664 
665  // vertices
666  dof_map[idx(0,0,0)] = 0;
667  Nodes.IntPoint(0).Set3(0., 0., 0.);
668  dof_map[idx(p,0,0)] = 1;
669  Nodes.IntPoint(1).Set3(1., 0., 0.);
670  dof_map[idx(0,p,0)] = 2;
671  Nodes.IntPoint(2).Set3(0., 1., 0.);
672  dof_map[idx(0,0,p)] = 3;
673  Nodes.IntPoint(3).Set3(0., 0., 1.);
674 
675  // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
676  int o = 4;
677  for (int i = 1; i < p; i++) // (0,1)
678  {
679  dof_map[idx(i,0,0)] = o;
680  Nodes.IntPoint(o++).Set3(double(i)/p, 0., 0.);
681  }
682  for (int i = 1; i < p; i++) // (0,2)
683  {
684  dof_map[idx(0,i,0)] = o;
685  Nodes.IntPoint(o++).Set3(0., double(i)/p, 0.);
686  }
687  for (int i = 1; i < p; i++) // (0,3)
688  {
689  dof_map[idx(0,0,i)] = o;
690  Nodes.IntPoint(o++).Set3(0., 0., double(i)/p);
691  }
692  for (int i = 1; i < p; i++) // (1,2)
693  {
694  dof_map[idx(p-i,i,0)] = o;
695  Nodes.IntPoint(o++).Set3(double(p-i)/p, double(i)/p, 0.);
696  }
697  for (int i = 1; i < p; i++) // (1,3)
698  {
699  dof_map[idx(p-i,0,i)] = o;
700  Nodes.IntPoint(o++).Set3(double(p-i)/p, 0., double(i)/p);
701  }
702  for (int i = 1; i < p; i++) // (2,3)
703  {
704  dof_map[idx(0,p-i,i)] = o;
705  Nodes.IntPoint(o++).Set3(0., double(p-i)/p, double(i)/p);
706  }
707 
708  // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
709  for (int j = 1; j < p; j++)
710  for (int i = 1; i + j < p; i++) // (1,2,3)
711  {
712  dof_map[idx(p-i-j,i,j)] = o;
713  Nodes.IntPoint(o++).Set3(double(p-i-j)/p, double(i)/p, double(j)/p);
714  }
715  for (int j = 1; j < p; j++)
716  for (int i = 1; i + j < p; i++) // (0,3,2)
717  {
718  dof_map[idx(0,j,i)] = o;
719  Nodes.IntPoint(o++).Set3(0., double(j)/p, double(i)/p);
720  }
721  for (int j = 1; j < p; j++)
722  for (int i = 1; i + j < p; i++) // (0,1,3)
723  {
724  dof_map[idx(i,0,j)] = o;
725  Nodes.IntPoint(o++).Set3(double(i)/p, 0., double(j)/p);
726  }
727  for (int j = 1; j < p; j++)
728  for (int i = 1; i + j < p; i++) // (0,2,1)
729  {
730  dof_map[idx(j,i,0)] = o;
731  Nodes.IntPoint(o++).Set3(double(j)/p, double(i)/p, 0.);
732  }
733 
734  // interior
735  for (int k = 1; k < p; k++)
736  for (int j = 1; j + k < p; j++)
737  for (int i = 1; i + j + k < p; i++)
738  {
739  dof_map[idx(i,j,k)] = o;
740  Nodes.IntPoint(o++).Set3(double(i)/p, double(j)/p, double(k)/p);
741  }
742 }
743 
744 // static method
746  const int p, const double l1, const double l2, const double l3,
747  double *shape)
748 {
749  const double l4 = 1. - l1 - l2 - l3;
750 
751  // The basis functions are the terms in the expansion:
752  // (l1 + l2 + l3 + l4)^p =
753  // \sum_{k=0}^p \binom{p}{k} l3^k
754  // \sum_{j=0}^{p-k} \binom{p-k}{j} l2^j
755  // \sum_{i=0}^{p-k-j} \binom{p-k-j}{i} l1^i l4^{p-k-j-i}
756  const int *bp = Poly_1D::Binom(p);
757  double l3k = 1.;
758  for (int o = 0, k = 0; k <= p; k++)
759  {
760  const int *bpk = Poly_1D::Binom(p - k);
761  const double ek = bp[k]*l3k;
762  double l2j = 1.;
763  for (int j = 0; j <= p - k; j++)
764  {
765  Poly_1D::CalcBinomTerms(p - k - j, l1, l4, &shape[o]);
766  double ekj = ek*bpk[j]*l2j;
767  for (int i = 0; i <= p - k - j; i++)
768  {
769  shape[o++] *= ekj;
770  }
771  l2j *= l2;
772  }
773  l3k *= l3;
774  }
775 }
776 
777 // static method
779  const int p, const double l1, const double l2, const double l3,
780  double *dshape_1d, double *dshape)
781 {
782  const int dof = ((p + 1)*(p + 2)*(p + 3))/6;
783  const double l4 = 1. - l1 - l2 - l3;
784 
785  // For the x derivatives, differentiate the terms of the expression:
786  // \sum_{k=0}^p \binom{p}{k} l3^k
787  // \sum_{j=0}^{p-k} \binom{p-k}{j} l2^j
788  // \sum_{i=0}^{p-k-j} \binom{p-k-j}{i} l1^i l4^{p-k-j-i}
789  const int *bp = Poly_1D::Binom(p);
790  double l3k = 1.;
791  for (int o = 0, k = 0; k <= p; k++)
792  {
793  const int *bpk = Poly_1D::Binom(p - k);
794  const double ek = bp[k]*l3k;
795  double l2j = 1.;
796  for (int j = 0; j <= p - k; j++)
797  {
798  Poly_1D::CalcDBinomTerms(p - k - j, l1, l4, dshape_1d);
799  double ekj = ek*bpk[j]*l2j;
800  for (int i = 0; i <= p - k - j; i++)
801  {
802  dshape[o++] = dshape_1d[i]*ekj;
803  }
804  l2j *= l2;
805  }
806  l3k *= l3;
807  }
808  // For the y derivatives, differentiate the terms of the expression:
809  // \sum_{k=0}^p \binom{p}{k} l3^k
810  // \sum_{i=0}^{p-k} \binom{p-k}{i} l1^i
811  // \sum_{j=0}^{p-k-i} \binom{p-k-i}{j} l2^j l4^{p-k-j-i}
812  l3k = 1.;
813  for (int ok = 0, k = 0; k <= p; k++)
814  {
815  const int *bpk = Poly_1D::Binom(p - k);
816  const double ek = bp[k]*l3k;
817  double l1i = 1.;
818  for (int i = 0; i <= p - k; i++)
819  {
820  Poly_1D::CalcDBinomTerms(p - k - i, l2, l4, dshape_1d);
821  double eki = ek*bpk[i]*l1i;
822  int o = ok + i;
823  for (int j = 0; j <= p - k - i; j++)
824  {
825  dshape[dof + o] = dshape_1d[j]*eki;
826  o += p - k - j + 1;
827  }
828  l1i *= l1;
829  }
830  l3k *= l3;
831  ok += ((p - k + 2)*(p - k + 1))/2;
832  }
833  // For the z derivatives, differentiate the terms of the expression:
834  // \sum_{j=0}^p \binom{p}{j} l2^j
835  // \sum_{i=0}^{p-j} \binom{p-j}{i} l1^i
836  // \sum_{k=0}^{p-j-i} \binom{p-j-i}{k} l3^k l4^{p-k-j-i}
837  double l2j = 1.;
838  for (int j = 0; j <= p; j++)
839  {
840  const int *bpj = Poly_1D::Binom(p - j);
841  const double ej = bp[j]*l2j;
842  double l1i = 1.;
843  for (int i = 0; i <= p - j; i++)
844  {
845  Poly_1D::CalcDBinomTerms(p - j - i, l3, l4, dshape_1d);
846  double eji = ej*bpj[i]*l1i;
847  int m = ((p + 2)*(p + 1))/2;
848  int n = ((p - j + 2)*(p - j + 1))/2;
849  for (int o = i, k = 0; k <= p - j - i; k++)
850  {
851  // m = ((p - k + 2)*(p - k + 1))/2;
852  // n = ((p - k - j + 2)*(p - k - j + 1))/2;
853  o += m;
854  dshape[2*dof + o - n] = dshape_1d[k]*eji;
855  m -= p - k + 1;
856  n -= p - k - j + 1;
857  }
858  l1i *= l1;
859  }
860  l2j *= l2;
861  }
862 }
863 
865  Vector &shape) const
866 {
867 #ifdef MFEM_THREAD_SAFE
868  Vector m_shape(dof);
869 #endif
870  CalcShape(order, ip.x, ip.y, ip.z, m_shape.GetData());
871  for (int i = 0; i < dof; i++)
872  {
873  shape(dof_map[i]) = m_shape(i);
874  }
875 }
876 
878  DenseMatrix &dshape) const
879 {
880 #ifdef MFEM_THREAD_SAFE
881  Vector dshape_1d(order + 1);
883 #endif
884  CalcDShape(order, ip.x, ip.y, ip.z, dshape_1d.GetData(), m_dshape.Data());
885  for (int d = 0; d < 3; d++)
886  {
887  for (int i = 0; i < dof; i++)
888  {
889  dshape(dof_map[i],d) = m_dshape(i,d);
890  }
891  }
892 }
893 
894 
896  : PositiveFiniteElement(3, Geometry::PRISM,
897  ((p + 1)*(p + 1)*(p + 2))/2, p, FunctionSpace::Qk),
898  TriangleFE(p),
899  SegmentFE(p)
900 {
901 #ifndef MFEM_THREAD_SAFE
906 #endif
907 
908  t_dof.SetSize(dof);
909  s_dof.SetSize(dof);
910 
911  // Nodal DoFs
912  t_dof[0] = 0; s_dof[0] = 0;
913  t_dof[1] = 1; s_dof[1] = 0;
914  t_dof[2] = 2; s_dof[2] = 0;
915  t_dof[3] = 0; s_dof[3] = 1;
916  t_dof[4] = 1; s_dof[4] = 1;
917  t_dof[5] = 2; s_dof[5] = 1;
918 
919  // Edge DoFs
920  int ne = p-1;
921  for (int i=1; i<p; i++)
922  {
923  t_dof[5 + 0 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 0 * ne + i] = 0;
924  t_dof[5 + 1 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 1 * ne + i] = 0;
925  t_dof[5 + 2 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 2 * ne + i] = 0;
926  t_dof[5 + 3 * ne + i] = 2 + 0 * ne + i; s_dof[5 + 3 * ne + i] = 1;
927  t_dof[5 + 4 * ne + i] = 2 + 1 * ne + i; s_dof[5 + 4 * ne + i] = 1;
928  t_dof[5 + 5 * ne + i] = 2 + 2 * ne + i; s_dof[5 + 5 * ne + i] = 1;
929  t_dof[5 + 6 * ne + i] = 0; s_dof[5 + 6 * ne + i] = i + 1;
930  t_dof[5 + 7 * ne + i] = 1; s_dof[5 + 7 * ne + i] = i + 1;
931  t_dof[5 + 8 * ne + i] = 2; s_dof[5 + 8 * ne + i] = i + 1;
932  }
933 
934  // Triangular Face DoFs
935  int k=0;
936  int nt = (p-1)*(p-2)/2;
937  for (int j=1; j<p; j++)
938  {
939  for (int i=1; i<j; i++)
940  {
941  t_dof[6 + 9 * ne + k] = 3 * p + k; s_dof[6 + 9 * ne + k] = 0;
942  t_dof[6 + 9 * ne + nt + k] = 3 * p + k; s_dof[6 + 9 * ne + nt + k] = 1;
943  k++;
944  }
945  }
946 
947  // Quadrilateral Face DoFs
948  k=0;
949  int nq = (p-1)*(p-1);
950  for (int j=1; j<p; j++)
951  {
952  for (int i=1; i<p; i++)
953  {
954  t_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 2 + 0 * ne + i;
955  t_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 2 + 1 * ne + i;
956  t_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 2 + 2 * ne + i;
957 
958  s_dof[6 + 9 * ne + 2 * nt + 0 * nq + k] = 1 + j;
959  s_dof[6 + 9 * ne + 2 * nt + 1 * nq + k] = 1 + j;
960  s_dof[6 + 9 * ne + 2 * nt + 2 * nq + k] = 1 + j;
961 
962  k++;
963  }
964  }
965 
966  // Interior DoFs
967  int m=0;
968  for (k=1; k<p; k++)
969  {
970  int l=0;
971  for (int j=1; j<p; j++)
972  {
973  for (int i=1; i<j; i++)
974  {
975  t_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 3 * p + l;
976  s_dof[6 + 9 * ne + 2 * nt + 3 * nq + m] = 1 + k;
977  l++; m++;
978  }
979  }
980  }
981 
982  // Define Nodes
983  const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
984  const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
985  for (int i=0; i<dof; i++)
986  {
987  Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
988  Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
989  Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
990  }
991 }
992 
994  Vector &shape) const
995 {
996 #ifdef MFEM_THREAD_SAFE
999 #endif
1000 
1001  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1002 
1004  SegmentFE.CalcShape(ipz, s_shape);
1005 
1006  for (int i=0; i<dof; i++)
1007  {
1008  shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
1009  }
1010 }
1011 
1013  DenseMatrix &dshape) const
1014 {
1015 #ifdef MFEM_THREAD_SAFE
1020 #endif
1021 
1022  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1023 
1026  SegmentFE.CalcShape(ipz, s_shape);
1028 
1029  for (int i=0; i<dof; i++)
1030  {
1031  dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
1032  dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
1033  dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
1034  }
1035 }
1036 
1038  : PositiveTensorFiniteElement(1, p, L2_DOF_MAP)
1039 {
1040 #ifndef MFEM_THREAD_SAFE
1041  shape_x.SetSize(p + 1);
1042  dshape_x.SetDataAndSize(NULL, p + 1);
1043 #endif
1044 
1045  if (p == 0)
1046  {
1047  Nodes.IntPoint(0).x = 0.5;
1048  }
1049  else
1050  {
1051  for (int i = 0; i <= p; i++)
1052  {
1053  Nodes.IntPoint(i).x = double(i)/p;
1054  }
1055  }
1056 }
1057 
1059  Vector &shape) const
1060 {
1061  Poly_1D::CalcBernstein(order, ip.x, shape);
1062 }
1063 
1065  DenseMatrix &dshape) const
1066 {
1067 #ifdef MFEM_THREAD_SAFE
1068  Vector shape_x(dof), dshape_x(dshape.Data(), dof);
1069 #else
1070  dshape_x.SetData(dshape.Data());
1071 #endif
1072  Poly_1D::CalcBernstein(order, ip.x, shape_x, dshape_x);
1073 }
1074 
1075 void L2Pos_SegmentElement::ProjectDelta(int vertex, Vector &dofs) const
1076 {
1077  dofs = 0.0;
1078  dofs[vertex*order] = 1.0;
1079 }
1080 
1081 
1083  : PositiveTensorFiniteElement(2, p, L2_DOF_MAP)
1084 {
1085 #ifndef MFEM_THREAD_SAFE
1086  shape_x.SetSize(p + 1);
1087  shape_y.SetSize(p + 1);
1088  dshape_x.SetSize(p + 1);
1089  dshape_y.SetSize(p + 1);
1090 #endif
1091 
1092  if (p == 0)
1093  {
1094  Nodes.IntPoint(0).Set2(0.5, 0.5);
1095  }
1096  else
1097  {
1098  for (int o = 0, j = 0; j <= p; j++)
1099  for (int i = 0; i <= p; i++)
1100  {
1101  Nodes.IntPoint(o++).Set2(double(i)/p, double(j)/p);
1102  }
1103  }
1104 }
1105 
1107  Vector &shape) const
1108 {
1109  const int p = order;
1110 
1111 #ifdef MFEM_THREAD_SAFE
1112  Vector shape_x(p+1), shape_y(p+1);
1113 #endif
1114 
1115  Poly_1D::CalcBernstein(p, ip.x, shape_x);
1116  Poly_1D::CalcBernstein(p, ip.y, shape_y);
1117 
1118  for (int o = 0, j = 0; j <= p; j++)
1119  for (int i = 0; i <= p; i++)
1120  {
1121  shape(o++) = shape_x(i)*shape_y(j);
1122  }
1123 }
1124 
1126  DenseMatrix &dshape) const
1127 {
1128  const int p = order;
1129 
1130 #ifdef MFEM_THREAD_SAFE
1131  Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
1132 #endif
1133 
1134  Poly_1D::CalcBernstein(p, ip.x, shape_x, dshape_x);
1135  Poly_1D::CalcBernstein(p, ip.y, shape_y, dshape_y);
1136 
1137  for (int o = 0, j = 0; j <= p; j++)
1138  for (int i = 0; i <= p; i++)
1139  {
1140  dshape(o,0) = dshape_x(i)* shape_y(j);
1141  dshape(o,1) = shape_x(i)*dshape_y(j); o++;
1142  }
1143 }
1144 
1146 {
1147  const int p = order;
1148 
1149  dofs = 0.0;
1150  switch (vertex)
1151  {
1152  case 0: dofs[0] = 1.0; break;
1153  case 1: dofs[p] = 1.0; break;
1154  case 2: dofs[p*(p + 2)] = 1.0; break;
1155  case 3: dofs[p*(p + 1)] = 1.0; break;
1156  }
1157 }
1158 
1159 
1161  : PositiveTensorFiniteElement(3, p, L2_DOF_MAP)
1162 {
1163 #ifndef MFEM_THREAD_SAFE
1164  shape_x.SetSize(p + 1);
1165  shape_y.SetSize(p + 1);
1166  shape_z.SetSize(p + 1);
1167  dshape_x.SetSize(p + 1);
1168  dshape_y.SetSize(p + 1);
1169  dshape_z.SetSize(p + 1);
1170 #endif
1171 
1172  if (p == 0)
1173  {
1174  Nodes.IntPoint(0).Set3(0.5, 0.5, 0.5);
1175  }
1176  else
1177  {
1178  for (int o = 0, k = 0; k <= p; k++)
1179  for (int j = 0; j <= p; j++)
1180  for (int i = 0; i <= p; i++)
1181  {
1182  Nodes.IntPoint(o++).Set3(double(i)/p, double(j)/p, double(k)/p);
1183  }
1184  }
1185 }
1186 
1188  Vector &shape) const
1189 {
1190  const int p = order;
1191 
1192 #ifdef MFEM_THREAD_SAFE
1193  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
1194 #endif
1195 
1196  Poly_1D::CalcBernstein(p, ip.x, shape_x);
1197  Poly_1D::CalcBernstein(p, ip.y, shape_y);
1198  Poly_1D::CalcBernstein(p, ip.z, shape_z);
1199 
1200  for (int o = 0, k = 0; k <= p; k++)
1201  for (int j = 0; j <= p; j++)
1202  for (int i = 0; i <= p; i++)
1203  {
1204  shape(o++) = shape_x(i)*shape_y(j)*shape_z(k);
1205  }
1206 }
1207 
1209  DenseMatrix &dshape) const
1210 {
1211  const int p = order;
1212 
1213 #ifdef MFEM_THREAD_SAFE
1214  Vector shape_x(p+1), shape_y(p+1), shape_z(p+1);
1215  Vector dshape_x(p+1), dshape_y(p+1), dshape_z(p+1);
1216 #endif
1217 
1218  Poly_1D::CalcBernstein(p, ip.x, shape_x, dshape_x);
1219  Poly_1D::CalcBernstein(p, ip.y, shape_y, dshape_y);
1220  Poly_1D::CalcBernstein(p, ip.z, shape_z, dshape_z);
1221 
1222  for (int o = 0, k = 0; k <= p; k++)
1223  for (int j = 0; j <= p; j++)
1224  for (int i = 0; i <= p; i++)
1225  {
1226  dshape(o,0) = dshape_x(i)* shape_y(j)* shape_z(k);
1227  dshape(o,1) = shape_x(i)*dshape_y(j)* shape_z(k);
1228  dshape(o,2) = shape_x(i)* shape_y(j)*dshape_z(k); o++;
1229  }
1230 }
1231 
1232 void L2Pos_HexahedronElement::ProjectDelta(int vertex, Vector &dofs) const
1233 {
1234  const int p = order;
1235 
1236  dofs = 0.0;
1237  switch (vertex)
1238  {
1239  case 0: dofs[0] = 1.0; break;
1240  case 1: dofs[p] = 1.0; break;
1241  case 2: dofs[p*(p + 2)] = 1.0; break;
1242  case 3: dofs[p*(p + 1)] = 1.0; break;
1243  case 4: dofs[p*(p + 1)*(p + 1)] = 1.0; break;
1244  case 5: dofs[p + p*(p + 1)*(p + 1)] = 1.0; break;
1245  case 6: dofs[dof - 1] = 1.0; break;
1246  case 7: dofs[dof - p - 1] = 1.0; break;
1247  }
1248 }
1249 
1250 
1252  : PositiveFiniteElement(2, Geometry::TRIANGLE, ((p + 1)*(p + 2))/2, p,
1253  FunctionSpace::Pk)
1254 {
1255 #ifndef MFEM_THREAD_SAFE
1256  dshape_1d.SetSize(p + 1);
1257 #endif
1258 
1259  if (p == 0)
1260  {
1261  Nodes.IntPoint(0).Set2(1./3, 1./3);
1262  }
1263  else
1264  {
1265  for (int o = 0, j = 0; j <= p; j++)
1266  for (int i = 0; i + j <= p; i++)
1267  {
1268  Nodes.IntPoint(o++).Set2(double(i)/p, double(j)/p);
1269  }
1270  }
1271 }
1272 
1274  Vector &shape) const
1275 {
1276  H1Pos_TriangleElement::CalcShape(order, ip.x, ip.y, shape.GetData());
1277 }
1278 
1280  DenseMatrix &dshape) const
1281 {
1282 #ifdef MFEM_THREAD_SAFE
1283  Vector dshape_1d(order + 1);
1284 #endif
1285 
1286  H1Pos_TriangleElement::CalcDShape(order, ip.x, ip.y, dshape_1d.GetData(),
1287  dshape.Data());
1288 }
1289 
1290 void L2Pos_TriangleElement::ProjectDelta(int vertex, Vector &dofs) const
1291 {
1292  dofs = 0.0;
1293  switch (vertex)
1294  {
1295  case 0: dofs[0] = 1.0; break;
1296  case 1: dofs[order] = 1.0; break;
1297  case 2: dofs[dof-1] = 1.0; break;
1298  }
1299 }
1300 
1301 
1303  : PositiveFiniteElement(3, Geometry::TETRAHEDRON,
1304  ((p + 1)*(p + 2)*(p + 3))/6, p, FunctionSpace::Pk)
1305 {
1306 #ifndef MFEM_THREAD_SAFE
1307  dshape_1d.SetSize(p + 1);
1308 #endif
1309 
1310  if (p == 0)
1311  {
1312  Nodes.IntPoint(0).Set3(0.25, 0.25, 0.25);
1313  }
1314  else
1315  {
1316  for (int o = 0, k = 0; k <= p; k++)
1317  for (int j = 0; j + k <= p; j++)
1318  for (int i = 0; i + j + k <= p; i++)
1319  {
1320  Nodes.IntPoint(o++).Set3(double(i)/p, double(j)/p, double(k)/p);
1321  }
1322  }
1323 }
1324 
1326  Vector &shape) const
1327 {
1329  shape.GetData());
1330 }
1331 
1333  DenseMatrix &dshape) const
1334 {
1335 #ifdef MFEM_THREAD_SAFE
1336  Vector dshape_1d(order + 1);
1337 #endif
1338 
1340  dshape_1d.GetData(), dshape.Data());
1341 }
1342 
1343 void L2Pos_TetrahedronElement::ProjectDelta(int vertex, Vector &dofs) const
1344 {
1345  dofs = 0.0;
1346  switch (vertex)
1347  {
1348  case 0: dofs[0] = 1.0; break;
1349  case 1: dofs[order] = 1.0; break;
1350  case 2: dofs[(order*(order+3))/2] = 1.0; break;
1351  case 3: dofs[dof-1] = 1.0; break;
1352  }
1353 }
1354 
1355 
1357  : PositiveFiniteElement(3, Geometry::PRISM,
1358  ((p + 1)*(p + 1)*(p + 2))/2, p, FunctionSpace::Qk),
1359  TriangleFE(p),
1360  SegmentFE(p)
1361 {
1362 #ifndef MFEM_THREAD_SAFE
1367 #endif
1368 
1369  t_dof.SetSize(dof);
1370  s_dof.SetSize(dof);
1371 
1372  // Interior DoFs
1373  int m=0;
1374  for (int k=0; k<=p; k++)
1375  {
1376  int l=0;
1377  for (int j=0; j<=p; j++)
1378  {
1379  for (int i=0; i<=j; i++)
1380  {
1381  t_dof[m] = l;
1382  s_dof[m] = k;
1383  l++; m++;
1384  }
1385  }
1386  }
1387 
1388  // Define Nodes
1389  const IntegrationRule & t_Nodes = TriangleFE.GetNodes();
1390  const IntegrationRule & s_Nodes = SegmentFE.GetNodes();
1391  for (int i=0; i<dof; i++)
1392  {
1393  Nodes.IntPoint(i).x = t_Nodes.IntPoint(t_dof[i]).x;
1394  Nodes.IntPoint(i).y = t_Nodes.IntPoint(t_dof[i]).y;
1395  Nodes.IntPoint(i).z = s_Nodes.IntPoint(s_dof[i]).x;
1396  }
1397 }
1398 
1400  Vector &shape) const
1401 {
1402 #ifdef MFEM_THREAD_SAFE
1405 #endif
1406 
1407  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1408 
1410  SegmentFE.CalcShape(ipz, s_shape);
1411 
1412  for (int i=0; i<dof; i++)
1413  {
1414  shape[i] = t_shape[t_dof[i]] * s_shape[s_dof[i]];
1415  }
1416 }
1417 
1419  DenseMatrix &dshape) const
1420 {
1421 #ifdef MFEM_THREAD_SAFE
1426 #endif
1427 
1428  IntegrationPoint ipz; ipz.x = ip.z; ipz.y = 0.0; ipz.z = 0.0;
1429 
1432  SegmentFE.CalcShape(ipz, s_shape);
1434 
1435  for (int i=0; i<dof; i++)
1436  {
1437  dshape(i, 0) = t_dshape(t_dof[i],0) * s_shape[s_dof[i]];
1438  dshape(i, 1) = t_dshape(t_dof[i],1) * s_shape[s_dof[i]];
1439  dshape(i, 2) = t_shape[t_dof[i]] * s_dshape(s_dof[i],0);
1440  }
1441 }
1442 
1443 }
Abstract class for all finite elements.
Definition: fe_base.hpp:235
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:349
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
L2Pos_WedgeElement(const int p)
Construct the L2Pos_WedgeElement of order p.
Definition: fe_pos.cpp:1356
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:1232
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:1145
L2Pos_TriangleElement(const int p)
Construct the L2Pos_TriangleElement of order p.
Definition: fe_pos.cpp:1251
H1Pos_SegmentElement(const int p)
Construct the H1Pos_SegmentElement of order p.
Definition: fe_pos.cpp:291
Base class for vector Coefficients that optionally depend on time and space.
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:1290
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 ...
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_pos.cpp:110
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_pos.cpp:329
PositiveTensorFiniteElement(const int dims, const int p, const DofMapType dmtype)
Definition: fe_pos.cpp:79
Class for finite elements utilizing the always positive Bernstein basis.
Definition: fe_pos.hpp:22
L2Pos_TetrahedronElement(const int p)
Construct the L2Pos_TetrahedronElement of order p.
Definition: fe_pos.cpp:1302
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_pos.cpp:1012
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_base.cpp:656
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_pos.cpp:134
H1Pos_SegmentElement SegmentFE
Definition: fe_pos.hpp:245
H1Pos_TetrahedronElement(const int p)
Construct the H1Pos_TetrahedronElement of order p.
Definition: fe_pos.cpp:643
int dim
Dimension of reference space.
Definition: fe_base.hpp:238
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
BiQuadPos2DFiniteElement()
Construct the BiQuadPos2DFiniteElement.
Definition: fe_pos.cpp:87
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
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_pos.cpp:1418
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_pos.cpp:1332
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_pos.cpp:1058
H1Pos_HexahedronElement(const int p)
Construct the H1Pos_HexahedronElement of order p.
Definition: fe_pos.cpp:424
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
static void CalcShape(const int p, const double x, const double y, const double z, double *shape)
Definition: fe_pos.cpp:745
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_pos.cpp:396
Possible basis types. Note that not all elements can use all BasisType(s).
Definition: fe_base.hpp:24
L2Pos_HexahedronElement(const int p)
Construct the L2Pos_HexahedronElement of order p.
Definition: fe_pos.cpp:1160
static void CalcDShape(const int p, const double x, const double y, const double z, double *dshape_1d, double *dshape)
Definition: fe_pos.cpp:778
static void CalcBinomTerms(const int p, const double x, const double y, double *u)
Compute the p terms in the expansion of the binomial (x + y)^p and store them in the already allocate...
Definition: fe_base.cpp:1938
H1Pos_TriangleElement(const int p)
Construct the H1Pos_TriangleElement of order p.
Definition: fe_pos.cpp:498
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition: fe_pos.cpp:183
Class for standard nodal finite elements.
Definition: fe_base.hpp:706
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_pos.cpp:280
static void CalcShape(const int p, const double x, const double y, double *shape)
Definition: fe_pos.cpp:553
Array< int > s_dof
Definition: fe_pos.hpp:356
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:417
H1Pos_WedgeElement(const int p)
Construct the H1Pos_WedgeElement of order p.
Definition: fe_pos.cpp:895
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_pos.cpp:1279
Array< int > t_dof
Definition: fe_pos.hpp:356
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:640
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
L2Pos_QuadrilateralElement(const int p)
Construct the L2Pos_QuadrilateralElement of order p.
Definition: fe_pos.cpp:1082
static void CalcDShape(const int p, const double x, const double y, double *dshape_1d, double *dshape)
Definition: fe_pos.cpp:579
void Set2(const double x1, const double x2)
Definition: intrules.hpp:80
int GetVDim()
Returns dimension of the vector.
void SetData(double *d)
Definition: vector.hpp:150
Array< int > t_dof
Definition: fe_pos.hpp:242
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_pos.cpp:1273
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_pos.cpp:1106
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
IntegrationRule Nodes
Definition: fe_base.hpp:248
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_pos.cpp:1187
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_pos.cpp:309
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:491
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
H1Pos_TriangleElement TriangleFE
Definition: fe_pos.hpp:244
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
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_pos.cpp:1325
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
DenseMatrix t_dshape
Definition: fe_pos.hpp:240
DenseMatrix t_dshape
Definition: fe_pos.hpp:354
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients &quot;p choose k&quot; for k=0,...,p for the given p.
Definition: fe_base.cpp:1889
DenseMatrix s_dshape
Definition: fe_pos.hpp:240
DenseMatrix s_dshape
Definition: fe_pos.hpp:354
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3924
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_pos.cpp:1399
L2Pos_TriangleElement TriangleFE
Definition: fe_pos.hpp:358
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:157
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
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_pos.cpp:446
H1Pos_QuadrilateralElement(const int p)
Construct the H1Pos_QuadrilateralElement of order p.
Definition: fe_pos.cpp:356
Class for integration point with weight.
Definition: intrules.hpp:25
L2Pos_SegmentElement SegmentFE
Definition: fe_pos.hpp:359
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_pos.cpp:467
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_pos.cpp:1064
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_pos.cpp:1208
int dof
Number of degrees of freedom.
Definition: fe_base.hpp:245
Array< int > s_dof
Definition: fe_pos.hpp:242
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:1075
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
L2Pos_SegmentElement(const int p)
Construct the L2Pos_SegmentElement of order p.
Definition: fe_pos.cpp:1037
Vector data type.
Definition: vector.hpp:60
static void CalcBernstein(const int p, const double x, double *u)
Compute the values of the Bernstein basis functions of order p at coordinate x and store the results ...
Definition: fe_base.hpp:1135
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
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_pos.cpp:376
Describes the function space on each element.
Definition: fe_base.hpp:218
RefCoord s[3]
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_pos.cpp:1125
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_pos.cpp:212
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
QuadPos1DFiniteElement()
Construct the QuadPos1DFiniteElement.
Definition: fe_pos.cpp:262
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_pos.cpp:993
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
static void CalcDBinomTerms(const int p, const double x, const double y, double *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p assuming t...
Definition: fe_base.cpp:2002
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_pos.cpp:270
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:245
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_pos.cpp:1343
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_pos.cpp:23