MFEM  v4.6.0
Finite element discretization library
lininteg.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #include "fem.hpp"
13 #include <cmath>
14 
15 namespace mfem
16 {
17 
19  const Array<int> &markers,
20  Vector &b)
21 {
22  MFEM_ABORT("Not supported.");
23 }
24 
26  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
27 {
28  mfem_error("LinearFormIntegrator::AssembleRHSElementVect(...)");
29 }
30 
32  const FiniteElement &el1, const FiniteElement &el2,
34 {
35  mfem_error("LinearFormIntegrator::AssembleRHSElementVect(...)");
36 }
37 
40  Vector &elvect)
41 {
42  int dof = el.GetDof();
43 
44  shape.SetSize(dof); // vector of size dof
45  elvect.SetSize(dof);
46  elvect = 0.0;
47 
48  const IntegrationRule *ir = IntRule;
49  if (ir == NULL)
50  {
51  // ir = &IntRules.Get(el.GetGeomType(),
52  // oa * el.GetOrder() + ob + Tr.OrderW());
53  ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
54  }
55 
56  for (int i = 0; i < ir->GetNPoints(); i++)
57  {
58  const IntegrationPoint &ip = ir->IntPoint(i);
59 
60  Tr.SetIntPoint (&ip);
61  double val = Tr.Weight() * Q.Eval(Tr, ip);
62 
63  el.CalcPhysShape(Tr, shape);
64 
65  add(elvect, ip.weight * val, shape, elvect);
66  }
67 }
68 
70  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
71 {
72  MFEM_ASSERT(delta != NULL, "coefficient must be DeltaCoefficient");
73  elvect.SetSize(fe.GetDof());
74  fe.CalcPhysShape(Trans, elvect);
75  elvect *= delta->EvalDelta(Trans, Trans.GetIntPoint());
76 }
77 
79  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
80 {
81  int dof = el.GetDof();
82  int spaceDim = Tr.GetSpaceDim();
83 
84  dshape.SetSize(dof, spaceDim);
85 
86  elvect.SetSize(dof);
87  elvect = 0.0;
88 
89  const IntegrationRule *ir = IntRule;
90  if (ir == NULL)
91  {
92  int intorder = 2 * el.GetOrder();
93  ir = &IntRules.Get(el.GetGeomType(), intorder);
94  }
95 
96  for (int i = 0; i < ir->GetNPoints(); i++)
97  {
98  const IntegrationPoint &ip = ir->IntPoint(i);
99 
100  Tr.SetIntPoint(&ip);
101  el.CalcPhysDShape(Tr, dshape);
102 
103  Q.Eval(Qvec, Tr, ip);
104  Qvec *= ip.weight * Tr.Weight();
105 
106  dshape.AddMult(Qvec, elvect);
107  }
108 }
109 
111  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
112 {
113  MFEM_ASSERT(vec_delta != NULL,"coefficient must be VectorDeltaCoefficient");
114  int dof = fe.GetDof();
115  int spaceDim = Trans.GetSpaceDim();
116 
117  dshape.SetSize(dof, spaceDim);
118  fe.CalcPhysDShape(Trans, dshape);
119 
120  vec_delta->EvalDelta(Qvec, Trans, Trans.GetIntPoint());
121 
122  elvect.SetSize(dof);
123  dshape.Mult(Qvec, elvect);
124 }
125 
127  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
128 {
129  int dof = el.GetDof();
130 
131  shape.SetSize(dof); // vector of size dof
132  elvect.SetSize(dof);
133  elvect = 0.0;
134 
135  const IntegrationRule *ir = IntRule;
136  if (ir == NULL)
137  {
138  int intorder = oa * el.GetOrder() + ob; // <----------
139  ir = &IntRules.Get(el.GetGeomType(), intorder);
140  }
141 
142  for (int i = 0; i < ir->GetNPoints(); i++)
143  {
144  const IntegrationPoint &ip = ir->IntPoint(i);
145 
146  Tr.SetIntPoint (&ip);
147  double val = Tr.Weight() * Q.Eval(Tr, ip);
148 
149  el.CalcShape(ip, shape);
150 
151  add(elvect, ip.weight * val, shape, elvect);
152  }
153 }
154 
156  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
157 {
158  int dof = el.GetDof();
159 
160  shape.SetSize(dof); // vector of size dof
161  elvect.SetSize(dof);
162  elvect = 0.0;
163 
164  const IntegrationRule *ir = IntRule;
165  if (ir == NULL)
166  {
167  int intorder = oa * el.GetOrder() + ob; // <------ user control
168  ir = &IntRules.Get(Tr.FaceGeom, intorder); // of integration order
169  }
170 
171  for (int i = 0; i < ir->GetNPoints(); i++)
172  {
173  const IntegrationPoint &ip = ir->IntPoint(i);
174 
175  // Set the integration point in the face and the neighboring element
176  Tr.SetAllIntPoints(&ip);
177 
178  // Access the neighboring element's integration point
179  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
180 
181  double val = Tr.Face->Weight() * ip.weight * Q.Eval(*Tr.Face, ip);
182 
183  el.CalcShape(eip, shape);
184 
185  add(elvect, val, shape, elvect);
186  }
187 }
188 
190  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
191 {
192  int dim = el.GetDim()+1;
193  int dof = el.GetDof();
194  Vector nor(dim), Qvec;
195 
196  shape.SetSize(dof);
197  elvect.SetSize(dof);
198  elvect = 0.0;
199 
200  const IntegrationRule *ir = IntRule;
201  if (ir == NULL)
202  {
203  int intorder = oa * el.GetOrder() + ob; // <----------
204  ir = &IntRules.Get(el.GetGeomType(), intorder);
205  }
206 
207  for (int i = 0; i < ir->GetNPoints(); i++)
208  {
209  const IntegrationPoint &ip = ir->IntPoint(i);
210 
211  Tr.SetIntPoint(&ip);
212  if (dim > 1)
213  {
214  CalcOrtho(Tr.Jacobian(), nor);
215  }
216  else
217  {
218  nor[0] = 1.0;
219  }
220  Q.Eval(Qvec, Tr, ip);
221 
222  el.CalcShape(ip, shape);
223 
224  elvect.Add(ip.weight*(Qvec*nor), shape);
225  }
226 }
227 
229  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
230 {
231  int dim = el.GetDim()+1;
232  int dof = el.GetDof();
233  Vector tangent(dim), Qvec;
234 
235  shape.SetSize(dof);
236  elvect.SetSize(dof);
237  elvect = 0.0;
238 
239  if (dim != 2)
240  {
241  mfem_error("These methods make sense only in 2D problems.");
242  }
243 
244  const IntegrationRule *ir = IntRule;
245  if (ir == NULL)
246  {
247  int intorder = oa * el.GetOrder() + ob; // <----------
248  ir = &IntRules.Get(el.GetGeomType(), intorder);
249  }
250 
251  for (int i = 0; i < ir->GetNPoints(); i++)
252  {
253  const IntegrationPoint &ip = ir->IntPoint(i);
254 
255  Tr.SetIntPoint(&ip);
256  const DenseMatrix &Jac = Tr.Jacobian();
257  tangent(0) = Jac(0,0);
258  tangent(1) = Jac(1,0);
259 
260  Q.Eval(Qvec, Tr, ip);
261 
262  el.CalcShape(ip, shape);
263 
264  add(elvect, ip.weight*(Qvec*tangent), shape, elvect);
265  }
266 }
267 
269  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
270 {
271  int vdim = Q.GetVDim();
272  int dof = el.GetDof();
273 
274  double val,cf;
275 
276  shape.SetSize(dof); // vector of size dof
277 
278  elvect.SetSize(dof * vdim);
279  elvect = 0.0;
280 
281  const IntegrationRule *ir = IntRule;
282  if (ir == NULL)
283  {
284  int intorder = 2*el.GetOrder();
285  ir = &IntRules.Get(el.GetGeomType(), intorder);
286  }
287 
288  for (int i = 0; i < ir->GetNPoints(); i++)
289  {
290  const IntegrationPoint &ip = ir->IntPoint(i);
291 
292  Tr.SetIntPoint (&ip);
293  val = Tr.Weight();
294 
295  el.CalcShape(ip, shape);
296  Q.Eval (Qvec, Tr, ip);
297 
298  for (int k = 0; k < vdim; k++)
299  {
300  cf = val * Qvec(k);
301 
302  for (int s = 0; s < dof; s++)
303  {
304  elvect(dof*k+s) += ip.weight * cf * shape(s);
305  }
306  }
307  }
308 }
309 
311  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
312 {
313  MFEM_ASSERT(vec_delta != NULL, "coefficient must be VectorDeltaCoefficient");
314  int vdim = Q.GetVDim();
315  int dof = fe.GetDof();
316 
317  shape.SetSize(dof);
318  fe.CalcPhysShape(Trans, shape);
319 
320  vec_delta->EvalDelta(Qvec, Trans, Trans.GetIntPoint());
321 
322  elvect.SetSize(dof*vdim);
323  DenseMatrix elvec_as_mat(elvect.GetData(), dof, vdim);
324  MultVWt(shape, Qvec, elvec_as_mat);
325 }
326 
328  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
329 {
330  const int dim = el.GetDim();
331  const int dof = el.GetDof();
332  const int vdim = Q.GetVDim();
333  const int sdim = Tr.GetSpaceDim();
334 
335  dshape.SetSize(dof,sdim);
336 
337  elvect.SetSize(dof*(vdim/sdim));
338  elvect = 0.0;
339 
340  const IntegrationRule *ir = IntRule;
341  if (ir == NULL)
342  {
343  int intorder = 2 * el.GetOrder();
344  ir = &IntRules.Get(el.GetGeomType(), intorder);
345  }
346 
347  Vector pelvect(dof);
348  Vector part_x(dim);
349 
350  for (int q = 0; q < ir->GetNPoints(); q++)
351  {
352  const IntegrationPoint &ip = ir->IntPoint(q);
353 
354  Tr.SetIntPoint(&ip);
355  el.CalcPhysDShape(Tr, dshape);
356 
357  Q.Eval(Qvec, Tr, ip);
358  Qvec *= ip.weight * Tr.Weight();
359 
360  for (int k = 0; k < vdim/sdim; k++)
361  {
362  for (int d=0; d < sdim; ++d) { part_x(d) = Qvec(k*sdim+d); }
363  dshape.Mult(part_x, pelvect);
364  for (int s = 0; s < dof; ++s) { elvect(s+k*dof) += pelvect(s); }
365  }
366  }
367 }
368 
371 {
372  MFEM_ABORT("Not implemented!");
373 }
374 
376  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
377 {
378  int vdim = Q.GetVDim();
379  int dof = el.GetDof();
380 
381  shape.SetSize(dof);
382  vec.SetSize(vdim);
383 
384  elvect.SetSize(dof * vdim);
385  elvect = 0.0;
386 
387  const IntegrationRule *ir = IntRule;
388  if (ir == NULL)
389  {
390  int intorder = 2*el.GetOrder();
391  ir = &IntRules.Get(el.GetGeomType(), intorder);
392  }
393 
394  for (int i = 0; i < ir->GetNPoints(); i++)
395  {
396  const IntegrationPoint &ip = ir->IntPoint(i);
397 
398  Q.Eval(vec, Tr, ip);
399  Tr.SetIntPoint (&ip);
400  vec *= Tr.Weight() * ip.weight;
401  el.CalcShape(ip, shape);
402  for (int k = 0; k < vdim; k++)
403  for (int s = 0; s < dof; s++)
404  {
405  elvect(dof*k+s) += vec(k) * shape(s);
406  }
407  }
408 }
409 
411  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
412 {
413  int vdim = Q.GetVDim();
414  int dof = el.GetDof();
415 
416  shape.SetSize(dof);
417  vec.SetSize(vdim);
418 
419  elvect.SetSize(dof * vdim);
420  elvect = 0.0;
421 
422  const IntegrationRule *ir = IntRule;
423  if (ir == NULL)
424  {
425  int intorder = 2*el.GetOrder();
426  ir = &IntRules.Get(Tr.GetGeometryType(), intorder);
427  }
428 
429  for (int i = 0; i < ir->GetNPoints(); i++)
430  {
431  const IntegrationPoint &ip = ir->IntPoint(i);
432 
433  // Set the integration point in the face and the neighboring element
434  Tr.SetAllIntPoints(&ip);
435 
436  // Access the neighboring element's integration point
437  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
438 
439  // Use Tr transformation in case Q depends on boundary attribute
440  Q.Eval(vec, Tr, ip);
441  vec *= Tr.Weight() * ip.weight;
442  el.CalcShape(eip, shape);
443  for (int k = 0; k < vdim; k++)
444  {
445  for (int s = 0; s < dof; s++)
446  {
447  elvect(dof*k+s) += vec(k) * shape(s);
448  }
449  }
450  }
451 }
452 
454  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
455 {
456  int dof = el.GetDof();
457  int spaceDim = Tr.GetSpaceDim();
458  int vdim = std::max(spaceDim, el.GetVDim());
459 
460  vshape.SetSize(dof,vdim);
461  vec.SetSize(vdim);
462 
463  elvect.SetSize(dof);
464  elvect = 0.0;
465 
466  const IntegrationRule *ir = IntRule;
467  if (ir == NULL)
468  {
469  // int intorder = 2*el.GetOrder() - 1; // ok for O(h^{k+1}) conv. in L2
470  int intorder = 2*el.GetOrder();
471  ir = &IntRules.Get(el.GetGeomType(), intorder);
472  }
473 
474  for (int i = 0; i < ir->GetNPoints(); i++)
475  {
476  const IntegrationPoint &ip = ir->IntPoint(i);
477 
478  Tr.SetIntPoint (&ip);
479  el.CalcVShape(Tr, vshape);
480 
481  QF.Eval (vec, Tr, ip);
482  vec *= ip.weight * Tr.Weight();
483  vshape.AddMult (vec, elvect);
484  }
485 }
486 
488  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
489 {
490  MFEM_ASSERT(vec_delta != NULL, "coefficient must be VectorDeltaCoefficient");
491  int dof = fe.GetDof();
492  int spaceDim = Trans.GetSpaceDim();
493 
494  vshape.SetSize(dof, spaceDim);
495  fe.CalcPhysVShape(Trans, vshape);
496 
497  vec_delta->EvalDelta(vec, Trans, Trans.GetIntPoint());
498 
499  elvect.SetSize(dof);
500  vshape.Mult(vec, elvect);
501 }
502 
504  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
505 {
506  int dof = el.GetDof();
507  int spaceDim = Tr.GetSpaceDim();
508  int n=(spaceDim == 3)? spaceDim : 1;
509  curlshape.SetSize(dof,n);
510  vec.SetSize(n);
511 
512  elvect.SetSize(dof);
513  elvect = 0.0;
514 
515  const IntegrationRule *ir = IntRule;
516  if (ir == NULL)
517  {
518  int intorder = 2*el.GetOrder();
519  ir = &IntRules.Get(el.GetGeomType(), intorder);
520  }
521 
522  for (int i = 0; i < ir->GetNPoints(); i++)
523  {
524  const IntegrationPoint &ip = ir->IntPoint(i);
525 
526  Tr.SetIntPoint (&ip);
527  el.CalcPhysCurlShape(Tr, curlshape);
528  QF->Eval(vec, Tr, ip);
529 
530  vec *= ip.weight * Tr.Weight();
531  curlshape.AddMult (vec, elvect);
532  }
533 }
534 
536  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
537 {
538  int spaceDim = Trans.GetSpaceDim();
539  MFEM_ASSERT(vec_delta != NULL,
540  "coefficient must be VectorDeltaCoefficient");
541  int dof = fe.GetDof();
542  int n=(spaceDim == 3)? spaceDim : 1;
543  vec.SetSize(n);
544  curlshape.SetSize(dof, n);
545  elvect.SetSize(dof);
546  fe.CalcPhysCurlShape(Trans, curlshape);
547 
548  vec_delta->EvalDelta(vec, Trans, Trans.GetIntPoint());
549  curlshape.Mult(vec, elvect);
550 }
551 
553  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
554 {
555  int dof = el.GetDof();
556 
557  divshape.SetSize(dof); // vector of size dof
558  elvect.SetSize(dof);
559  elvect = 0.0;
560 
561  const IntegrationRule *ir = IntRule;
562  if (ir == NULL)
563  {
564  int intorder = 2 * el.GetOrder();
565  ir = &IntRules.Get(el.GetGeomType(), intorder);
566  }
567 
568  for (int i = 0; i < ir->GetNPoints(); i++)
569  {
570  const IntegrationPoint &ip = ir->IntPoint(i);
571 
572  Tr.SetIntPoint (&ip);
573  double val = Tr.Weight() * Q.Eval(Tr, ip);
574  el.CalcPhysDivShape(Tr, divshape);
575 
576  add(elvect, ip.weight * val, divshape, elvect);
577  }
578 }
579 
581  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
582 {
583  MFEM_ASSERT(delta != NULL, "coefficient must be DeltaCoefficient");
584  elvect.SetSize(fe.GetDof());
585  fe.CalcPhysDivShape(Trans, elvect);
586  elvect *= delta->EvalDelta(Trans, Trans.GetIntPoint());
587 }
588 
590  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
591 {
592  int dim = el.GetDim()+1;
593  int dof = el.GetDof();
594 
595  shape.SetSize (dof);
596  nor.SetSize (dim);
597  elvect.SetSize (dim*dof);
598 
599  const IntegrationRule *ir = IntRule;
600  if (ir == NULL)
601  {
602  ir = &IntRules.Get(el.GetGeomType(), el.GetOrder() + 1);
603  }
604 
605  elvect = 0.0;
606  for (int i = 0; i < ir->GetNPoints(); i++)
607  {
608  const IntegrationPoint &ip = ir->IntPoint(i);
609  Tr.SetIntPoint (&ip);
610  CalcOrtho(Tr.Jacobian(), nor);
611  el.CalcShape (ip, shape);
612  nor *= Sign * ip.weight * F -> Eval (Tr, ip);
613  for (int j = 0; j < dof; j++)
614  for (int k = 0; k < dim; k++)
615  {
616  elvect(dof*k+j) += nor(k) * shape(j);
617  }
618  }
619 }
620 
621 
623  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
624 {
625  int dof = el.GetDof();
626 
627  shape.SetSize(dof);
628  elvect.SetSize(dof);
629  elvect = 0.0;
630 
631  const IntegrationRule *ir = IntRule;
632  if (ir == NULL)
633  {
634  int intorder = oa * el.GetOrder() + ob; // <----------
635  ir = &IntRules.Get(el.GetGeomType(), intorder);
636  }
637 
638  for (int i = 0; i < ir->GetNPoints(); i++)
639  {
640  const IntegrationPoint &ip = ir->IntPoint(i);
641  el.CalcShape(ip, shape);
642 
643  double val = ip.weight;
644  if (F)
645  {
646  Tr.SetIntPoint (&ip);
647  val *= F->Eval(Tr, ip);
648  }
649 
650  elvect.Add(val, shape);
651  }
652 }
653 
655  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
656 {
657  int dof = el.GetDof();
658  int dim = el.GetDim();
659  int vdim = el.GetVDim();
660  DenseMatrix vshape(dof, vdim);
661  Vector f_loc(3);
662  Vector f_hat(2);
663 
664  MFEM_VERIFY(vdim == 2, "VectorFEBoundaryTangentLFIntegrator "
665  "must be called with vector basis functions of dimension 2.");
666 
667  elvect.SetSize(dof);
668  elvect = 0.0;
669 
670  const IntegrationRule *ir = IntRule;
671  if (ir == NULL)
672  {
673  int intorder = oa * el.GetOrder() + ob; // <----------
674  ir = &IntRules.Get(el.GetGeomType(), intorder);
675  }
676 
677  for (int i = 0; i < ir->GetNPoints(); i++)
678  {
679  const IntegrationPoint &ip = ir->IntPoint(i);
680 
681  el.CalcVShape(ip, vshape);
682 
683  Tr.SetIntPoint(&ip);
684  f.Eval(f_loc, Tr, ip);
685 
686  if (dim == 2)
687  {
688  Tr.Jacobian().MultTranspose(f_loc, f_hat);
689  }
690  else if (dim == 1)
691  {
692  const DenseMatrix & J = Tr.Jacobian();
693  f_hat(0) = J(0,0) * f_loc(0) + J(1,0) * f_loc(1);
694  f_hat(1) = f_loc(2);
695  }
696  else
697  {
698  f_hat(0) = f_loc(1);
699  f_hat(1) = f_loc(2);
700  }
701 
702  Swap<double>(f_hat(0), f_hat(1));
703  f_hat(0) = -f_hat(0);
704  f_hat *= ip.weight;
705 
706  vshape.AddMult(f_hat, elvect);
707  }
708 }
709 
711  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
712 {
713  mfem_error("BoundaryFlowIntegrator::AssembleRHSElementVect\n"
714  " is not implemented as boundary integrator!\n"
715  " Use LinearForm::AddBdrFaceIntegrator instead of\n"
716  " LinearForm::AddBoundaryIntegrator.");
717 }
718 
720  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
721 {
722  int dim, ndof, order;
723  double un, w, vu_data[3], nor_data[3];
724 
725  dim = el.GetDim();
726  ndof = el.GetDof();
727  Vector vu(vu_data, dim), nor(nor_data, dim);
728 
729  const IntegrationRule *ir = IntRule;
730  if (ir == NULL)
731  {
732  // Assuming order(u)==order(mesh)
733  order = Tr.Elem1->OrderW() + 2*el.GetOrder();
734  if (el.Space() == FunctionSpace::Pk)
735  {
736  order++;
737  }
738  ir = &IntRules.Get(Tr.GetGeometryType(), order);
739  }
740 
741  shape.SetSize(ndof);
742  elvect.SetSize(ndof);
743  elvect = 0.0;
744 
745  for (int p = 0; p < ir->GetNPoints(); p++)
746  {
747  const IntegrationPoint &ip = ir->IntPoint(p);
748 
749  // Set the integration point in the face and the neighboring element
750  Tr.SetAllIntPoints(&ip);
751 
752  // Access the neighboring element's integration point
753  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
754  el.CalcShape(eip, shape);
755 
756  // Use Tr.Elem1 transformation for u so that it matches the coefficient
757  // used with the ConvectionIntegrator and/or the DGTraceIntegrator.
758  u->Eval(vu, *Tr.Elem1, eip);
759 
760  if (dim == 1)
761  {
762  nor(0) = 2*eip.x - 1.0;
763  }
764  else
765  {
766  CalcOrtho(Tr.Jacobian(), nor);
767  }
768 
769  un = vu * nor;
770  w = 0.5*alpha*un - beta*fabs(un);
771  w *= ip.weight*f->Eval(Tr, ip);
772  elvect.Add(w, shape);
773  }
774 }
775 
777  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
778 {
779  mfem_error("DGDirichletLFIntegrator::AssembleRHSElementVect");
780 }
781 
783  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
784 {
785  int dim, ndof;
786  bool kappa_is_nonzero = (kappa != 0.);
787  double w;
788 
789  dim = el.GetDim();
790  ndof = el.GetDof();
791 
792  nor.SetSize(dim);
793  nh.SetSize(dim);
794  ni.SetSize(dim);
795  adjJ.SetSize(dim);
796  if (MQ)
797  {
798  mq.SetSize(dim);
799  }
800 
801  shape.SetSize(ndof);
802  dshape.SetSize(ndof, dim);
803  dshape_dn.SetSize(ndof);
804 
805  elvect.SetSize(ndof);
806  elvect = 0.0;
807 
808  const IntegrationRule *ir = IntRule;
809  if (ir == NULL)
810  {
811  // a simple choice for the integration order; is this OK?
812  int order = 2*el.GetOrder();
813  ir = &IntRules.Get(Tr.GetGeometryType(), order);
814  }
815 
816  for (int p = 0; p < ir->GetNPoints(); p++)
817  {
818  const IntegrationPoint &ip = ir->IntPoint(p);
819 
820  // Set the integration point in the face and the neighboring element
821  Tr.SetAllIntPoints(&ip);
822 
823  // Access the neighboring element's integration point
824  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
825 
826  if (dim == 1)
827  {
828  nor(0) = 2*eip.x - 1.0;
829  }
830  else
831  {
832  CalcOrtho(Tr.Jacobian(), nor);
833  }
834 
835  el.CalcShape(eip, shape);
836  el.CalcDShape(eip, dshape);
837 
838  // compute uD through the face transformation
839  w = ip.weight * uD->Eval(Tr, ip) / Tr.Elem1->Weight();
840  if (!MQ)
841  {
842  if (Q)
843  {
844  w *= Q->Eval(*Tr.Elem1, eip);
845  }
846  ni.Set(w, nor);
847  }
848  else
849  {
850  nh.Set(w, nor);
851  MQ->Eval(mq, *Tr.Elem1, eip);
852  mq.MultTranspose(nh, ni);
853  }
855  adjJ.Mult(ni, nh);
856 
858  elvect.Add(sigma, dshape_dn);
859 
860  if (kappa_is_nonzero)
861  {
862  elvect.Add(kappa*(ni*nor), shape);
863  }
864  }
865 }
866 
868  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
869 {
870  mfem_error("DGElasticityDirichletLFIntegrator::AssembleRHSElementVect");
871 }
872 
874  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
875 {
876  MFEM_ASSERT(Tr.Elem2No < 0, "interior boundary is not supported");
877 
878 #ifdef MFEM_THREAD_SAFE
879  Vector shape;
883  Vector nor;
886  Vector u_dir;
887 #endif
888 
889  const int dim = el.GetDim();
890  const int ndofs = el.GetDof();
891  const int nvdofs = dim*ndofs;
892 
893  elvect.SetSize(nvdofs);
894  elvect = 0.0;
895 
896  adjJ.SetSize(dim);
897  shape.SetSize(ndofs);
898  dshape.SetSize(ndofs, dim);
899  dshape_ps.SetSize(ndofs, dim);
900  nor.SetSize(dim);
901  dshape_dn.SetSize(ndofs);
902  dshape_du.SetSize(ndofs);
903  u_dir.SetSize(dim);
904 
905  const IntegrationRule *ir = IntRule;
906  if (ir == NULL)
907  {
908  const int order = 2*el.GetOrder(); // <-----
909  ir = &IntRules.Get(Tr.GetGeometryType(), order);
910  }
911 
912  for (int pi = 0; pi < ir->GetNPoints(); ++pi)
913  {
914  const IntegrationPoint &ip = ir->IntPoint(pi);
915 
916  // Set the integration point in the face and the neighboring element
917  Tr.SetAllIntPoints(&ip);
918 
919  // Access the neighboring element's integration point
920  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
921 
922  // Evaluate the Dirichlet b.c. using the face transformation.
923  uD.Eval(u_dir, Tr, ip);
924 
925  el.CalcShape(eip, shape);
926  el.CalcDShape(eip, dshape);
927 
930 
931  if (dim == 1)
932  {
933  nor(0) = 2*eip.x - 1.0;
934  }
935  else
936  {
937  CalcOrtho(Tr.Jacobian(), nor);
938  }
939 
940  double wL, wM, jcoef;
941  {
942  const double w = ip.weight / Tr.Elem1->Weight();
943  wL = w * lambda->Eval(*Tr.Elem1, eip);
944  wM = w * mu->Eval(*Tr.Elem1, eip);
945  jcoef = kappa * (wL + 2.0*wM) * (nor*nor);
948  }
949 
950  // alpha < uD, (lambda div(v) I + mu (grad(v) + grad(v)^T)) . n > +
951  // + kappa < h^{-1} (lambda + 2 mu) uD, v >
952 
953  // i = idof + ndofs * im
954  // v_phi(i,d) = delta(im,d) phi(idof)
955  // div(v_phi(i)) = dphi(idof,im)
956  // (grad(v_phi(i)))(k,l) = delta(im,k) dphi(idof,l)
957  //
958  // term 1:
959  // alpha < uD, lambda div(v_phi(i)) n >
960  // alpha lambda div(v_phi(i)) (uD.n) =
961  // alpha lambda dphi(idof,im) (uD.n) --> quadrature -->
962  // ip.weight/det(J1) alpha lambda (uD.nor) dshape_ps(idof,im) =
963  // alpha * wL * (u_dir*nor) * dshape_ps(idof,im)
964  // term 2:
965  // < alpha uD, mu grad(v_phi(i)).n > =
966  // alpha mu uD^T grad(v_phi(i)) n =
967  // alpha mu uD(k) delta(im,k) dphi(idof,l) n(l) =
968  // alpha mu uD(im) dphi(idof,l) n(l) --> quadrature -->
969  // ip.weight/det(J1) alpha mu uD(im) dshape_ps(idof,l) nor(l) =
970  // alpha * wM * u_dir(im) * dshape_dn(idof)
971  // term 3:
972  // < alpha uD, mu (grad(v_phi(i)))^T n > =
973  // alpha mu n^T grad(v_phi(i)) uD =
974  // alpha mu n(k) delta(im,k) dphi(idof,l) uD(l) =
975  // alpha mu n(im) dphi(idof,l) uD(l) --> quadrature -->
976  // ip.weight/det(J1) alpha mu nor(im) dshape_ps(idof,l) uD(l) =
977  // alpha * wM * nor(im) * dshape_du(idof)
978  // term j:
979  // < kappa h^{-1} (lambda + 2 mu) uD, v_phi(i) > =
980  // kappa/h (lambda + 2 mu) uD(k) v_phi(i,k) =
981  // kappa/h (lambda + 2 mu) uD(k) delta(im,k) phi(idof) =
982  // kappa/h (lambda + 2 mu) uD(im) phi(idof) --> quadrature -->
983  // [ 1/h = |nor|/det(J1) ]
984  // ip.weight/det(J1) |nor|^2 kappa (lambda + 2 mu) uD(im) phi(idof) =
985  // jcoef * u_dir(im) * shape(idof)
986 
987  wM *= alpha;
988  const double t1 = alpha * wL * (u_dir*nor);
989  for (int im = 0, i = 0; im < dim; ++im)
990  {
991  const double t2 = wM * u_dir(im);
992  const double t3 = wM * nor(im);
993  const double tj = jcoef * u_dir(im);
994  for (int idof = 0; idof < ndofs; ++idof, ++i)
995  {
996  elvect(i) += (t1*dshape_ps(idof,im) + t2*dshape_dn(idof) +
997  t3*dshape_du(idof) + tj*shape(idof));
998  }
999  }
1000  }
1001 }
1002 
1003 
1004 
1006 (const FiniteElement &el,
1008  Vector &elvect)
1009 {
1010  int n = el.GetDof();
1011  elvect.SetSize(n);
1012  for (int i = 0; i < n; i++)
1013  {
1014  elvect(i) = dist(generator);
1015  }
1016 
1017  int iel = Tr.ElementNo;
1018 
1019  if (!save_factors || !L[iel])
1020  {
1021  DenseMatrix *M, m;
1022  if (save_factors)
1023  {
1024  L[iel]=new DenseMatrix;
1025  M = L[iel];
1026  }
1027  else
1028  {
1029  M = &m;
1030  }
1031  massinteg.AssembleElementMatrix(el, Tr, *M);
1032  CholeskyFactors chol(M->Data());
1033  chol.Factor(M->Height());
1034  chol.LMult(n,1,elvect.GetData());
1035  }
1036  else
1037  {
1038  CholeskyFactors chol(L[iel]->Data());
1039  chol.LMult(n,1,elvect.GetData());
1040  }
1041 }
1042 
1043 
1045  const FiniteElement &fe, ElementTransformation &Tr, Vector &elvect)
1046 {
1047  const IntegrationRule *ir =
1049 
1050  const int nqp = ir->GetNPoints();
1051  const int vdim = vqfc.GetVDim();
1052  const int ndofs = fe.GetDof();
1053  Vector shape(ndofs);
1054  Vector temp(vdim);
1055  elvect.SetSize(vdim * ndofs);
1056  elvect = 0.0;
1057  for (int q = 0; q < nqp; q++)
1058  {
1059  const IntegrationPoint &ip = ir->IntPoint(q);
1060  Tr.SetIntPoint(&ip);
1061  const double w = Tr.Weight() * ip.weight;
1062  vqfc.Eval(temp, Tr, ip);
1063  fe.CalcShape(ip, shape);
1064  for (int ind = 0; ind < vdim; ind++)
1065  {
1066  for (int nd = 0; nd < ndofs; nd++)
1067  {
1068  elvect(nd + ind * ndofs) += w * shape(nd) * temp(ind);
1069  }
1070  }
1071  }
1072 }
1073 
1074 
1077  Vector &elvect)
1078 {
1079  const IntegrationRule *ir =
1081 
1082  const int nqp = ir->GetNPoints();
1083  const int ndofs = fe.GetDof();
1084  Vector shape(ndofs);
1085  elvect.SetSize(ndofs);
1086  elvect = 0.0;
1087  for (int q = 0; q < nqp; q++)
1088  {
1089  const IntegrationPoint &ip = ir->IntPoint(q);
1090  Tr.SetIntPoint (&ip);
1091  const double w = Tr.Weight() * ip.weight;
1092  double temp = qfc.Eval(Tr, ip);
1093  fe.CalcShape(ip, shape);
1094  shape *= (w * temp);
1095  elvect += shape;
1096  }
1097 }
1098 
1099 }
Abstract class for all finite elements.
Definition: fe_base.hpp:233
virtual void AssembleDeltaElementVect(const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
Assemble the delta coefficient at the IntegrationPoint set in Trans which is assumed to map to the de...
Definition: lininteg.cpp:110
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:622
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
virtual void AssembleDeltaElementVect(const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
Assemble the delta coefficient at the IntegrationPoint set in Trans which is assumed to map to the de...
Definition: lininteg.cpp:580
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:1006
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:3127
virtual int OrderW() const =0
Return the order of the determinant of the Jacobian (weight) of the transformation.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
const QuadratureFunction & GetQuadFunction() const
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:552
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:517
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:710
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe_base.hpp:337
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:867
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:78
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
virtual void AssembleDeltaElementVect(const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
Assemble the delta coefficient at the IntegrationPoint set in Trans which is assumed to map to the de...
Definition: lininteg.cpp:69
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2478
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:268
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
void CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape) const
Equivalent to the CalcVShape() method with the same arguments.
Definition: fe_base.hpp:411
const IntegrationRule * IntRule
Definition: lininteg.hpp:27
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition: eltrans.hpp:98
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2673
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:126
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:317
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
const QuadratureFunction & GetQuadFunction() const
Polynomials of order k.
Definition: fe_base.hpp:221
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
virtual void AssembleDeltaElementVect(const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect) override
Assemble the delta coefficient at the IntegrationPoint set in Trans which is assumed to map to the de...
Definition: lininteg.cpp:369
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:162
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:214
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
double b
Definition: lissajous.cpp:42
void SetAllIntPoints(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.hpp:566
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:375
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...
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_base.cpp:71
VectorDeltaCoefficient * vec_delta
Definition: lininteg.hpp:66
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
virtual void AssembleDeltaElementVect(const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
Assemble the delta coefficient at the IntegrationPoint set in Trans which is assumed to map to the de...
Definition: lininteg.cpp:535
MatrixCoefficient * MQ
Definition: lininteg.hpp:540
virtual void AssembleDevice(const FiniteElementSpace &fes, const Array< int > &markers, Vector &b)
Method defining assembly on device.
Definition: lininteg.cpp:18
int GetVDim()
Returns dimension of the vector.
void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in physical space at the po...
Definition: fe_base.cpp:58
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:172
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition: fe_base.cpp:192
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
void LMult(int m, int n, double *X) const
Definition: densemat.cpp:3646
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:228
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:776
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)=0
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
y += a * A.x
Definition: densemat.cpp:251
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
virtual bool Factor(int m, double TOL=0.0)
Compute the Cholesky factorization of the current matrix.
Definition: densemat.cpp:3594
DeltaCoefficient * delta
Definition: lininteg.hpp:65
const IntegrationPoint & GetElement1IntPoint()
Get a const reference to the integration point in neighboring element 1 corresponding to the currentl...
Definition: eltrans.hpp:577
const IntegrationRule & GetIntRule(int idx) const
Return the IntegrationRule associated with entity idx.
Definition: qspace.hpp:73
virtual void EvalDelta(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Return the specified direction vector multiplied by the value returned by DeltaCoefficient::EvalDelta...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:589
virtual double EvalDelta(ElementTransformation &T, const IntegrationPoint &ip)
The value of the function assuming we are evaluating at the delta center.
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:263
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:654
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition: fe_base.cpp:182
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:248
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:82
Class for integration point with weight.
Definition: intrules.hpp:31
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
Definition: fe_base.hpp:314
virtual void AssembleDeltaElementVect(const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
Assemble the delta coefficient at the IntegrationPoint set in Trans which is assumed to map to the de...
Definition: lininteg.cpp:487
int dim
Definition: ex24.cpp:53
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:40
virtual void AssembleRHSElementVect(const FiniteElement &fe, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:1044
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:189
ElementTransformation * Elem1
Definition: eltrans.hpp:522
ElementTransformation * Face
Definition: eltrans.hpp:523
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:503
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
virtual void AssembleRHSElementVect(const FiniteElement &fe, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:1075
Vector data type.
Definition: vector.hpp:58
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:453
RefCoord s[3]
virtual void AssembleDeltaElementVect(const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
Assemble the delta coefficient at the IntegrationPoint set in Trans which is assumed to map to the de...
Definition: lininteg.cpp:310
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:38
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect) override
Definition: lininteg.cpp:327
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 GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327