MFEM  v4.3.0 Finite element discretization library
lininteg.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12
13 #include "fem.hpp"
14 #include <cmath>
15
16 namespace mfem
17 {
18
20  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
21 {
22  mfem_error("LinearFormIntegrator::AssembleRHSElementVect(...)");
23 }
24
26  const FiniteElement &el1, const FiniteElement &el2,
28 {
29  mfem_error("LinearFormIntegrator::AssembleRHSElementVect(...)");
30 }
31
34  Vector &elvect)
35 {
36  int dof = el.GetDof();
37
38  shape.SetSize(dof); // vector of size dof
39  elvect.SetSize(dof);
40  elvect = 0.0;
41
42  const IntegrationRule *ir = IntRule;
43  if (ir == NULL)
44  {
45  // ir = &IntRules.Get(el.GetGeomType(),
46  // oa * el.GetOrder() + ob + Tr.OrderW());
47  ir = &IntRules.Get(el.GetGeomType(), oa * el.GetOrder() + ob);
48  }
49
50  for (int i = 0; i < ir->GetNPoints(); i++)
51  {
52  const IntegrationPoint &ip = ir->IntPoint(i);
53
54  Tr.SetIntPoint (&ip);
55  double val = Tr.Weight() * Q.Eval(Tr, ip);
56
57  el.CalcPhysShape(Tr, shape);
58
59  add(elvect, ip.weight * val, shape, elvect);
60  }
61 }
62
64  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
65 {
66  MFEM_ASSERT(delta != NULL, "coefficient must be DeltaCoefficient");
67  elvect.SetSize(fe.GetDof());
68  fe.CalcPhysShape(Trans, elvect);
69  elvect *= delta->EvalDelta(Trans, Trans.GetIntPoint());
70 }
71
73  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
74 {
75  int dof = el.GetDof();
76  int spaceDim = Tr.GetSpaceDim();
77
78  dshape.SetSize(dof, spaceDim);
79
80  elvect.SetSize(dof);
81  elvect = 0.0;
82
83  const IntegrationRule *ir = IntRule;
84  if (ir == NULL)
85  {
86  int intorder = 2 * el.GetOrder();
87  ir = &IntRules.Get(el.GetGeomType(), intorder);
88  }
89
90  for (int i = 0; i < ir->GetNPoints(); i++)
91  {
92  const IntegrationPoint &ip = ir->IntPoint(i);
93
94  Tr.SetIntPoint(&ip);
95  el.CalcPhysDShape(Tr, dshape);
96
97  Q.Eval(Qvec, Tr, ip);
98  Qvec *= ip.weight * Tr.Weight();
99
101  }
102 }
103
105  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
106 {
107  MFEM_ASSERT(vec_delta != NULL,"coefficient must be VectorDeltaCoefficient");
108  int dof = fe.GetDof();
109  int spaceDim = Trans.GetSpaceDim();
110
111  dshape.SetSize(dof, spaceDim);
112  fe.CalcPhysDShape(Trans, dshape);
113
114  vec_delta->EvalDelta(Qvec, Trans, Trans.GetIntPoint());
115
116  elvect.SetSize(dof);
117  dshape.Mult(Qvec, elvect);
118 }
119
121  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
122 {
123  int dof = el.GetDof();
124
125  shape.SetSize(dof); // vector of size dof
126  elvect.SetSize(dof);
127  elvect = 0.0;
128
129  const IntegrationRule *ir = IntRule;
130  if (ir == NULL)
131  {
132  int intorder = oa * el.GetOrder() + ob; // <----------
133  ir = &IntRules.Get(el.GetGeomType(), intorder);
134  }
135
136  for (int i = 0; i < ir->GetNPoints(); i++)
137  {
138  const IntegrationPoint &ip = ir->IntPoint(i);
139
140  Tr.SetIntPoint (&ip);
141  double val = Tr.Weight() * Q.Eval(Tr, ip);
142
143  el.CalcShape(ip, shape);
144
145  add(elvect, ip.weight * val, shape, elvect);
146  }
147 }
148
150  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
151 {
152  int dof = el.GetDof();
153
154  shape.SetSize(dof); // vector of size dof
155  elvect.SetSize(dof);
156  elvect = 0.0;
157
158  const IntegrationRule *ir = IntRule;
159  if (ir == NULL)
160  {
161  int intorder = oa * el.GetOrder() + ob; // <------ user control
162  ir = &IntRules.Get(Tr.FaceGeom, intorder); // of integration order
163  }
164
165  for (int i = 0; i < ir->GetNPoints(); i++)
166  {
167  const IntegrationPoint &ip = ir->IntPoint(i);
168
169  // Set the integration point in the face and the neighboring element
170  Tr.SetAllIntPoints(&ip);
171
172  // Access the neighboring element's integration point
173  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
174
175  double val = Tr.Face->Weight() * ip.weight * Q.Eval(*Tr.Face, ip);
176
177  el.CalcShape(eip, shape);
178
180  }
181 }
182
184  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
185 {
186  int dim = el.GetDim()+1;
187  int dof = el.GetDof();
188  Vector nor(dim), Qvec;
189
190  shape.SetSize(dof);
191  elvect.SetSize(dof);
192  elvect = 0.0;
193
194  const IntegrationRule *ir = IntRule;
195  if (ir == NULL)
196  {
197  int intorder = oa * el.GetOrder() + ob; // <----------
198  ir = &IntRules.Get(el.GetGeomType(), intorder);
199  }
200
201  for (int i = 0; i < ir->GetNPoints(); i++)
202  {
203  const IntegrationPoint &ip = ir->IntPoint(i);
204
205  Tr.SetIntPoint(&ip);
206  CalcOrtho(Tr.Jacobian(), nor);
207  Q.Eval(Qvec, Tr, ip);
208
209  el.CalcShape(ip, shape);
210
212  }
213 }
214
216  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
217 {
218  int dim = el.GetDim()+1;
219  int dof = el.GetDof();
220  Vector tangent(dim), Qvec;
221
222  shape.SetSize(dof);
223  elvect.SetSize(dof);
224  elvect = 0.0;
225
226  if (dim != 2)
227  {
228  mfem_error("These methods make sense only in 2D problems.");
229  }
230
231  const IntegrationRule *ir = IntRule;
232  if (ir == NULL)
233  {
234  int intorder = oa * el.GetOrder() + ob; // <----------
235  ir = &IntRules.Get(el.GetGeomType(), intorder);
236  }
237
238  for (int i = 0; i < ir->GetNPoints(); i++)
239  {
240  const IntegrationPoint &ip = ir->IntPoint(i);
241
242  Tr.SetIntPoint(&ip);
243  const DenseMatrix &Jac = Tr.Jacobian();
244  tangent(0) = Jac(0,0);
245  tangent(1) = Jac(1,0);
246
247  Q.Eval(Qvec, Tr, ip);
248
249  el.CalcShape(ip, shape);
250
252  }
253 }
254
256  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
257 {
258  int vdim = Q.GetVDim();
259  int dof = el.GetDof();
260
261  double val,cf;
262
263  shape.SetSize(dof); // vector of size dof
264
265  elvect.SetSize(dof * vdim);
266  elvect = 0.0;
267
268  const IntegrationRule *ir = IntRule;
269  if (ir == NULL)
270  {
271  int intorder = 2*el.GetOrder();
272  ir = &IntRules.Get(el.GetGeomType(), intorder);
273  }
274
275  for (int i = 0; i < ir->GetNPoints(); i++)
276  {
277  const IntegrationPoint &ip = ir->IntPoint(i);
278
279  Tr.SetIntPoint (&ip);
280  val = Tr.Weight();
281
282  el.CalcShape(ip, shape);
283  Q.Eval (Qvec, Tr, ip);
284
285  for (int k = 0; k < vdim; k++)
286  {
287  cf = val * Qvec(k);
288
289  for (int s = 0; s < dof; s++)
290  {
291  elvect(dof*k+s) += ip.weight * cf * shape(s);
292  }
293  }
294  }
295 }
296
298  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
299 {
300  MFEM_ASSERT(vec_delta != NULL, "coefficient must be VectorDeltaCoefficient");
301  int vdim = Q.GetVDim();
302  int dof = fe.GetDof();
303
304  shape.SetSize(dof);
305  fe.CalcPhysShape(Trans, shape);
306
307  vec_delta->EvalDelta(Qvec, Trans, Trans.GetIntPoint());
308
309  elvect.SetSize(dof*vdim);
310  DenseMatrix elvec_as_mat(elvect.GetData(), dof, vdim);
311  MultVWt(shape, Qvec, elvec_as_mat);
312 }
313
315  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
316 {
317  int vdim = Q.GetVDim();
318  int dof = el.GetDof();
319
320  shape.SetSize(dof);
321  vec.SetSize(vdim);
322
323  elvect.SetSize(dof * vdim);
324  elvect = 0.0;
325
326  const IntegrationRule *ir = IntRule;
327  if (ir == NULL)
328  {
329  int intorder = 2*el.GetOrder();
330  ir = &IntRules.Get(el.GetGeomType(), intorder);
331  }
332
333  for (int i = 0; i < ir->GetNPoints(); i++)
334  {
335  const IntegrationPoint &ip = ir->IntPoint(i);
336
337  Q.Eval(vec, Tr, ip);
338  Tr.SetIntPoint (&ip);
339  vec *= Tr.Weight() * ip.weight;
340  el.CalcShape(ip, shape);
341  for (int k = 0; k < vdim; k++)
342  for (int s = 0; s < dof; s++)
343  {
344  elvect(dof*k+s) += vec(k) * shape(s);
345  }
346  }
347 }
348
350  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
351 {
352  int vdim = Q.GetVDim();
353  int dof = el.GetDof();
354
355  shape.SetSize(dof);
356  vec.SetSize(vdim);
357
358  elvect.SetSize(dof * vdim);
359  elvect = 0.0;
360
361  const IntegrationRule *ir = IntRule;
362  if (ir == NULL)
363  {
364  int intorder = 2*el.GetOrder();
365  ir = &IntRules.Get(Tr.GetGeometryType(), intorder);
366  }
367
368  for (int i = 0; i < ir->GetNPoints(); i++)
369  {
370  const IntegrationPoint &ip = ir->IntPoint(i);
371
372  // Set the integration point in the face and the neighboring element
373  Tr.SetAllIntPoints(&ip);
374
375  // Access the neighboring element's integration point
376  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
377
378  // Use Tr transformation in case Q depends on boundary attribute
379  Q.Eval(vec, Tr, ip);
380  vec *= Tr.Weight() * ip.weight;
381  el.CalcShape(eip, shape);
382  for (int k = 0; k < vdim; k++)
383  {
384  for (int s = 0; s < dof; s++)
385  {
386  elvect(dof*k+s) += vec(k) * shape(s);
387  }
388  }
389  }
390 }
391
393  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
394 {
395  int dof = el.GetDof();
396  int spaceDim = Tr.GetSpaceDim();
397
398  vshape.SetSize(dof,spaceDim);
399  vec.SetSize(spaceDim);
400
401  elvect.SetSize(dof);
402  elvect = 0.0;
403
404  const IntegrationRule *ir = IntRule;
405  if (ir == NULL)
406  {
407  // int intorder = 2*el.GetOrder() - 1; // ok for O(h^{k+1}) conv. in L2
408  int intorder = 2*el.GetOrder();
409  ir = &IntRules.Get(el.GetGeomType(), intorder);
410  }
411
412  for (int i = 0; i < ir->GetNPoints(); i++)
413  {
414  const IntegrationPoint &ip = ir->IntPoint(i);
415
416  Tr.SetIntPoint (&ip);
417  el.CalcVShape(Tr, vshape);
418
419  QF.Eval (vec, Tr, ip);
420  vec *= ip.weight * Tr.Weight();
422  }
423 }
424
426  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
427 {
428  MFEM_ASSERT(vec_delta != NULL, "coefficient must be VectorDeltaCoefficient");
429  int dof = fe.GetDof();
430  int spaceDim = Trans.GetSpaceDim();
431
432  vshape.SetSize(dof, spaceDim);
433  fe.CalcPhysVShape(Trans, vshape);
434
435  vec_delta->EvalDelta(vec, Trans, Trans.GetIntPoint());
436
437  elvect.SetSize(dof);
438  vshape.Mult(vec, elvect);
439 }
440
442  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
443 {
444  int dof = el.GetDof();
445  int spaceDim = Tr.GetSpaceDim();
446  int n=(spaceDim == 3)? spaceDim : 1;
447  curlshape.SetSize(dof,n);
448  vec.SetSize(n);
449
450  elvect.SetSize(dof);
451  elvect = 0.0;
452
453  const IntegrationRule *ir = IntRule;
454  if (ir == NULL)
455  {
456  int intorder = 2*el.GetOrder();
457  ir = &IntRules.Get(el.GetGeomType(), intorder);
458  }
459
460  for (int i = 0; i < ir->GetNPoints(); i++)
461  {
462  const IntegrationPoint &ip = ir->IntPoint(i);
463
464  Tr.SetIntPoint (&ip);
465  el.CalcPhysCurlShape(Tr, curlshape);
466  QF->Eval(vec, Tr, ip);
467
468  vec *= ip.weight * Tr.Weight();
470  }
471 }
472
474  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
475 {
476  int spaceDim = Trans.GetSpaceDim();
477  MFEM_ASSERT(vec_delta != NULL,
478  "coefficient must be VectorDeltaCoefficient");
479  int dof = fe.GetDof();
480  int n=(spaceDim == 3)? spaceDim : 1;
481  vec.SetSize(n);
482  curlshape.SetSize(dof, n);
483  elvect.SetSize(dof);
484  fe.CalcPhysCurlShape(Trans, curlshape);
485
486  vec_delta->EvalDelta(vec, Trans, Trans.GetIntPoint());
487  curlshape.Mult(vec, elvect);
488 }
489
491  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
492 {
493  int dof = el.GetDof();
494
495  divshape.SetSize(dof); // vector of size dof
496  elvect.SetSize(dof);
497  elvect = 0.0;
498
499  const IntegrationRule *ir = IntRule;
500  if (ir == NULL)
501  {
502  int intorder = 2 * el.GetOrder();
503  ir = &IntRules.Get(el.GetGeomType(), intorder);
504  }
505
506  for (int i = 0; i < ir->GetNPoints(); i++)
507  {
508  const IntegrationPoint &ip = ir->IntPoint(i);
509
510  Tr.SetIntPoint (&ip);
511  double val = Tr.Weight() * Q.Eval(Tr, ip);
512  el.CalcPhysDivShape(Tr, divshape);
513
514  add(elvect, ip.weight * val, divshape, elvect);
515  }
516 }
517
519  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
520 {
521  MFEM_ASSERT(delta != NULL, "coefficient must be DeltaCoefficient");
522  elvect.SetSize(fe.GetDof());
523  fe.CalcPhysDivShape(Trans, elvect);
524  elvect *= delta->EvalDelta(Trans, Trans.GetIntPoint());
525 }
526
528  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
529 {
530  int dim = el.GetDim()+1;
531  int dof = el.GetDof();
532
533  shape.SetSize (dof);
534  nor.SetSize (dim);
535  elvect.SetSize (dim*dof);
536
537  const IntegrationRule *ir = IntRule;
538  if (ir == NULL)
539  {
540  ir = &IntRules.Get(el.GetGeomType(), el.GetOrder() + 1);
541  }
542
543  elvect = 0.0;
544  for (int i = 0; i < ir->GetNPoints(); i++)
545  {
546  const IntegrationPoint &ip = ir->IntPoint(i);
547  Tr.SetIntPoint (&ip);
548  CalcOrtho(Tr.Jacobian(), nor);
549  el.CalcShape (ip, shape);
550  nor *= Sign * ip.weight * F -> Eval (Tr, ip);
551  for (int j = 0; j < dof; j++)
552  for (int k = 0; k < dim; k++)
553  {
554  elvect(dof*k+j) += nor(k) * shape(j);
555  }
556  }
557 }
558
559
561  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
562 {
563  int dof = el.GetDof();
564
565  shape.SetSize(dof);
566  elvect.SetSize(dof);
567  elvect = 0.0;
568
569  const IntegrationRule *ir = IntRule;
570  if (ir == NULL)
571  {
572  int intorder = oa * el.GetOrder() + ob; // <----------
573  ir = &IntRules.Get(el.GetGeomType(), intorder);
574  }
575
576  for (int i = 0; i < ir->GetNPoints(); i++)
577  {
578  const IntegrationPoint &ip = ir->IntPoint(i);
579  el.CalcShape(ip, shape);
580
581  double val = ip.weight;
582  if (F)
583  {
584  Tr.SetIntPoint (&ip);
585  val *= F->Eval(Tr, ip);
586  }
587
589  }
590 }
591
593  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
594 {
595  int dof = el.GetDof();
596  DenseMatrix vshape(dof, 2);
597  Vector f_loc(3);
598  Vector f_hat(2);
599
600  elvect.SetSize(dof);
601  elvect = 0.0;
602
603  const IntegrationRule *ir = IntRule;
604  if (ir == NULL)
605  {
606  int intorder = oa * el.GetOrder() + ob; // <----------
607  ir = &IntRules.Get(el.GetGeomType(), intorder);
608  }
609
610  for (int i = 0; i < ir->GetNPoints(); i++)
611  {
612  const IntegrationPoint &ip = ir->IntPoint(i);
613
614  Tr.SetIntPoint(&ip);
615  f.Eval(f_loc, Tr, ip);
616  Tr.Jacobian().MultTranspose(f_loc, f_hat);
617  el.CalcVShape(ip, vshape);
618
619  Swap<double>(f_hat(0), f_hat(1));
620  f_hat(0) = -f_hat(0);
621  f_hat *= ip.weight;
623  }
624 }
625
627  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
628 {
629  mfem_error("BoundaryFlowIntegrator::AssembleRHSElementVect\n"
630  " is not implemented as boundary integrator!\n"
633 }
634
636  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
637 {
638  int dim, ndof, order;
639  double un, w, vu_data[3], nor_data[3];
640
641  dim = el.GetDim();
642  ndof = el.GetDof();
643  Vector vu(vu_data, dim), nor(nor_data, dim);
644
645  const IntegrationRule *ir = IntRule;
646  if (ir == NULL)
647  {
648  // Assuming order(u)==order(mesh)
649  order = Tr.Elem1->OrderW() + 2*el.GetOrder();
650  if (el.Space() == FunctionSpace::Pk)
651  {
652  order++;
653  }
654  ir = &IntRules.Get(Tr.GetGeometryType(), order);
655  }
656
657  shape.SetSize(ndof);
658  elvect.SetSize(ndof);
659  elvect = 0.0;
660
661  for (int p = 0; p < ir->GetNPoints(); p++)
662  {
663  const IntegrationPoint &ip = ir->IntPoint(p);
664
665  // Set the integration point in the face and the neighboring element
666  Tr.SetAllIntPoints(&ip);
667
668  // Access the neighboring element's integration point
669  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
670  el.CalcShape(eip, shape);
671
672  // Use Tr.Elem1 transformation for u so that it matches the coefficient
673  // used with the ConvectionIntegrator and/or the DGTraceIntegrator.
674  u->Eval(vu, *Tr.Elem1, eip);
675
676  if (dim == 1)
677  {
678  nor(0) = 2*eip.x - 1.0;
679  }
680  else
681  {
682  CalcOrtho(Tr.Jacobian(), nor);
683  }
684
685  un = vu * nor;
686  w = 0.5*alpha*un - beta*fabs(un);
687  w *= ip.weight*f->Eval(Tr, ip);
689  }
690 }
691
693  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
694 {
695  mfem_error("DGDirichletLFIntegrator::AssembleRHSElementVect");
696 }
697
699  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
700 {
701  int dim, ndof;
702  bool kappa_is_nonzero = (kappa != 0.);
703  double w;
704
705  dim = el.GetDim();
706  ndof = el.GetDof();
707
708  nor.SetSize(dim);
709  nh.SetSize(dim);
710  ni.SetSize(dim);
712  if (MQ)
713  {
714  mq.SetSize(dim);
715  }
716
717  shape.SetSize(ndof);
718  dshape.SetSize(ndof, dim);
719  dshape_dn.SetSize(ndof);
720
721  elvect.SetSize(ndof);
722  elvect = 0.0;
723
724  const IntegrationRule *ir = IntRule;
725  if (ir == NULL)
726  {
727  // a simple choice for the integration order; is this OK?
728  int order = 2*el.GetOrder();
729  ir = &IntRules.Get(Tr.GetGeometryType(), order);
730  }
731
732  for (int p = 0; p < ir->GetNPoints(); p++)
733  {
734  const IntegrationPoint &ip = ir->IntPoint(p);
735
736  // Set the integration point in the face and the neighboring element
737  Tr.SetAllIntPoints(&ip);
738
739  // Access the neighboring element's integration point
740  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
741
742  if (dim == 1)
743  {
744  nor(0) = 2*eip.x - 1.0;
745  }
746  else
747  {
748  CalcOrtho(Tr.Jacobian(), nor);
749  }
750
751  el.CalcShape(eip, shape);
752  el.CalcDShape(eip, dshape);
753
754  // compute uD through the face transformation
755  w = ip.weight * uD->Eval(Tr, ip) / Tr.Elem1->Weight();
756  if (!MQ)
757  {
758  if (Q)
759  {
760  w *= Q->Eval(*Tr.Elem1, eip);
761  }
762  ni.Set(w, nor);
763  }
764  else
765  {
766  nh.Set(w, nor);
767  MQ->Eval(mq, *Tr.Elem1, eip);
768  mq.MultTranspose(nh, ni);
769  }
772
775
776  if (kappa_is_nonzero)
777  {
779  }
780  }
781 }
782
784  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
785 {
786  mfem_error("DGElasticityDirichletLFIntegrator::AssembleRHSElementVect");
787 }
788
790  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
791 {
792  MFEM_ASSERT(Tr.Elem2No < 0, "interior boundary is not supported");
793
795  Vector shape;
799  Vector nor;
802  Vector u_dir;
803 #endif
804
805  const int dim = el.GetDim();
806  const int ndofs = el.GetDof();
807  const int nvdofs = dim*ndofs;
808
809  elvect.SetSize(nvdofs);
810  elvect = 0.0;
811
813  shape.SetSize(ndofs);
814  dshape.SetSize(ndofs, dim);
815  dshape_ps.SetSize(ndofs, dim);
816  nor.SetSize(dim);
817  dshape_dn.SetSize(ndofs);
818  dshape_du.SetSize(ndofs);
819  u_dir.SetSize(dim);
820
821  const IntegrationRule *ir = IntRule;
822  if (ir == NULL)
823  {
824  const int order = 2*el.GetOrder(); // <-----
825  ir = &IntRules.Get(Tr.GetGeometryType(), order);
826  }
827
828  for (int pi = 0; pi < ir->GetNPoints(); ++pi)
829  {
830  const IntegrationPoint &ip = ir->IntPoint(pi);
831
832  // Set the integration point in the face and the neighboring element
833  Tr.SetAllIntPoints(&ip);
834
835  // Access the neighboring element's integration point
836  const IntegrationPoint &eip = Tr.GetElement1IntPoint();
837
838  // Evaluate the Dirichlet b.c. using the face transformation.
839  uD.Eval(u_dir, Tr, ip);
840
841  el.CalcShape(eip, shape);
842  el.CalcDShape(eip, dshape);
843
846
847  if (dim == 1)
848  {
849  nor(0) = 2*eip.x - 1.0;
850  }
851  else
852  {
853  CalcOrtho(Tr.Jacobian(), nor);
854  }
855
856  double wL, wM, jcoef;
857  {
858  const double w = ip.weight / Tr.Elem1->Weight();
859  wL = w * lambda->Eval(*Tr.Elem1, eip);
860  wM = w * mu->Eval(*Tr.Elem1, eip);
861  jcoef = kappa * (wL + 2.0*wM) * (nor*nor);
862  dshape_ps.Mult(nor, dshape_dn);
863  dshape_ps.Mult(u_dir, dshape_du);
864  }
865
866  // alpha < uD, (lambda div(v) I + mu (grad(v) + grad(v)^T)) . n > +
867  // + kappa < h^{-1} (lambda + 2 mu) uD, v >
868
869  // i = idof + ndofs * im
870  // v_phi(i,d) = delta(im,d) phi(idof)
871  // div(v_phi(i)) = dphi(idof,im)
872  // (grad(v_phi(i)))(k,l) = delta(im,k) dphi(idof,l)
873  //
874  // term 1:
875  // alpha < uD, lambda div(v_phi(i)) n >
876  // alpha lambda div(v_phi(i)) (uD.n) =
877  // alpha lambda dphi(idof,im) (uD.n) --> quadrature -->
878  // ip.weight/det(J1) alpha lambda (uD.nor) dshape_ps(idof,im) =
879  // alpha * wL * (u_dir*nor) * dshape_ps(idof,im)
880  // term 2:
881  // < alpha uD, mu grad(v_phi(i)).n > =
882  // alpha mu uD^T grad(v_phi(i)) n =
883  // alpha mu uD(k) delta(im,k) dphi(idof,l) n(l) =
884  // alpha mu uD(im) dphi(idof,l) n(l) --> quadrature -->
885  // ip.weight/det(J1) alpha mu uD(im) dshape_ps(idof,l) nor(l) =
886  // alpha * wM * u_dir(im) * dshape_dn(idof)
887  // term 3:
888  // < alpha uD, mu (grad(v_phi(i)))^T n > =
889  // alpha mu n^T grad(v_phi(i)) uD =
890  // alpha mu n(k) delta(im,k) dphi(idof,l) uD(l) =
891  // alpha mu n(im) dphi(idof,l) uD(l) --> quadrature -->
892  // ip.weight/det(J1) alpha mu nor(im) dshape_ps(idof,l) uD(l) =
893  // alpha * wM * nor(im) * dshape_du(idof)
894  // term j:
895  // < kappa h^{-1} (lambda + 2 mu) uD, v_phi(i) > =
896  // kappa/h (lambda + 2 mu) uD(k) v_phi(i,k) =
897  // kappa/h (lambda + 2 mu) uD(k) delta(im,k) phi(idof) =
898  // kappa/h (lambda + 2 mu) uD(im) phi(idof) --> quadrature -->
899  // [ 1/h = |nor|/det(J1) ]
900  // ip.weight/det(J1) |nor|^2 kappa (lambda + 2 mu) uD(im) phi(idof) =
901  // jcoef * u_dir(im) * shape(idof)
902
903  wM *= alpha;
904  const double t1 = alpha * wL * (u_dir*nor);
905  for (int im = 0, i = 0; im < dim; ++im)
906  {
907  const double t2 = wM * u_dir(im);
908  const double t3 = wM * nor(im);
909  const double tj = jcoef * u_dir(im);
910  for (int idof = 0; idof < ndofs; ++idof, ++i)
911  {
912  elvect(i) += (t1*dshape_ps(idof,im) + t2*dshape_dn(idof) +
913  t3*dshape_du(idof) + tj*shape(idof));
914  }
915  }
916  }
917 }
918
920  const FiniteElement &fe, ElementTransformation &Tr, Vector &elvect)
921 {
922  const IntegrationRule *ir =
924
925  const int nqp = ir->GetNPoints();
926  const int vdim = vqfc.GetVDim();
927  const int ndofs = fe.GetDof();
928  Vector shape(ndofs);
929  Vector temp(vdim);
930  elvect.SetSize(vdim * ndofs);
931  elvect = 0.0;
932  for (int q = 0; q < nqp; q++)
933  {
934  const IntegrationPoint &ip = ir->IntPoint(q);
935  Tr.SetIntPoint(&ip);
936  const double w = Tr.Weight() * ip.weight;
937  vqfc.Eval(temp, Tr, ip);
938  fe.CalcShape(ip, shape);
939  for (int ind = 0; ind < vdim; ind++)
940  {
941  for (int nd = 0; nd < ndofs; nd++)
942  {
943  elvect(nd + ind * ndofs) += w * shape(nd) * temp(ind);
944  }
945  }
946  }
947 }
948
951  Vector &elvect)
952 {
953  const IntegrationRule *ir =
955
956  const int nqp = ir->GetNPoints();
957  const int ndofs = fe.GetDof();
958  Vector shape(ndofs);
959  elvect.SetSize(ndofs);
960  elvect = 0.0;
961  for (int q = 0; q < nqp; q++)
962  {
963  const IntegrationPoint &ip = ir->IntPoint(q);
964  Tr.SetIntPoint (&ip);
965  const double w = Tr.Weight() * ip.weight;
966  double temp = qfc.Eval(Tr, ip);
967  fe.CalcShape(ip, shape);
968  shape *= (w * temp);
969  elvect += shape;
970  }
971 }
972
973 }
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.hpp:243
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:104
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:560
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:317
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
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:518
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.cpp:40
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:2757
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:920
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:149
const IntegrationRule & GetElementIntRule(int idx) const
Get the IntegrationRule associated with mesh element idx.
Definition: fespace.hpp:950
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:490
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.cpp:58
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
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:626
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:783
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:72
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
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:63
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
Definition: densemat.cpp:2089
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:327
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:255
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe.hpp:337
const IntegrationRule * IntRule
Definition: lininteg.hpp:25
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition: eltrans.hpp:90
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:199
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2303
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:120
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:291
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.cpp:71
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:320
Polynomials of order k.
Definition: fe.hpp:226
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:203
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
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:153
void SetAllIntPoints(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.hpp:553
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:111
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:314
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...
Definition: gridfunc.hpp:775
VectorDeltaCoefficient * vec_delta
Definition: lininteg.hpp:55
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
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:473
MatrixCoefficient * MQ
Definition: lininteg.hpp:448
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.cpp:192
int GetVDim()
Returns dimension of the vector.
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 AddMult(const Vector &x, Vector &y) const
y += A.x
Definition: densemat.cpp:226
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:215
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:692
void CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape) const
Equivalent to the CalcVShape() method with the same arguments.
Definition: fe.hpp:412
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)=0
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:323
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.cpp:182
DeltaCoefficient * delta
Definition: lininteg.hpp:54
const IntegrationPoint & GetElement1IntPoint()
Get a const reference to the integration point in neighboring element 1 corresponding to the currentl...
Definition: eltrans.hpp:564
virtual void EvalDelta(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Return the specified direction vector multiplied by the value returned by DeltaCoefficient::EvalDelta...
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:527
virtual double EvalDelta(ElementTransformation &T, const IntegrationPoint &ip)
The value of the function assuming we are evaluating at the delta center.
Definition: coefficient.cpp:83
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:258
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:592
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:243
Class for integration point with weight.
Definition: intrules.hpp:25
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:425
int dim
Definition: ex24.cpp:53
virtual void AssembleRHSElementVect(const FiniteElement &fe, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:919
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:183
ElementTransformation * Elem1
Definition: eltrans.hpp:509
ElementTransformation * Face
Definition: eltrans.hpp:510
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:441
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:949
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:175
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:392
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:297
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:32
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...