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