MFEM  v3.4
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
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 
66 
68  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
69 {
70  int dof = el.GetDof();
71 
72  shape.SetSize(dof); // vector of size dof
73  elvect.SetSize(dof);
74  elvect = 0.0;
75 
76  const IntegrationRule *ir = IntRule;
77  if (ir == NULL)
78  {
79  int intorder = oa * el.GetOrder() + ob; // <----------
80  ir = &IntRules.Get(el.GetGeomType(), intorder);
81  }
82 
83  for (int i = 0; i < ir->GetNPoints(); i++)
84  {
85  const IntegrationPoint &ip = ir->IntPoint(i);
86 
87  Tr.SetIntPoint (&ip);
88  double val = Tr.Weight() * Q.Eval(Tr, ip);
89 
90  el.CalcShape(ip, shape);
91 
92  add(elvect, ip.weight * val, shape, elvect);
93  }
94 }
95 
97  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
98 {
99  int dim = el.GetDim()+1;
100  int dof = el.GetDof();
101  Vector nor(dim), Qvec;
102 
103  shape.SetSize(dof);
104  elvect.SetSize(dof);
105  elvect = 0.0;
106 
107  const IntegrationRule *ir = IntRule;
108  if (ir == NULL)
109  {
110  int intorder = oa * el.GetOrder() + ob; // <----------
111  ir = &IntRules.Get(el.GetGeomType(), intorder);
112  }
113 
114  for (int i = 0; i < ir->GetNPoints(); i++)
115  {
116  const IntegrationPoint &ip = ir->IntPoint(i);
117 
118  Tr.SetIntPoint(&ip);
119  CalcOrtho(Tr.Jacobian(), nor);
120  Q.Eval(Qvec, Tr, ip);
121 
122  el.CalcShape(ip, shape);
123 
124  elvect.Add(ip.weight*(Qvec*nor), shape);
125  }
126 }
127 
129  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
130 {
131  int dim = el.GetDim()+1;
132  int dof = el.GetDof();
133  Vector tangent(dim), Qvec;
134 
135  shape.SetSize(dof);
136  elvect.SetSize(dof);
137  elvect = 0.0;
138 
139  if (dim != 2)
140  {
141  mfem_error("These methods make sense only in 2D problems.");
142  }
143 
144  const IntegrationRule *ir = IntRule;
145  if (ir == NULL)
146  {
147  int intorder = oa * el.GetOrder() + ob; // <----------
148  ir = &IntRules.Get(el.GetGeomType(), intorder);
149  }
150 
151  for (int i = 0; i < ir->GetNPoints(); i++)
152  {
153  const IntegrationPoint &ip = ir->IntPoint(i);
154 
155  Tr.SetIntPoint(&ip);
156  const DenseMatrix &Jac = Tr.Jacobian();
157  tangent(0) = Jac(0,0);
158  tangent(1) = Jac(1,0);
159 
160  Q.Eval(Qvec, Tr, ip);
161 
162  el.CalcShape(ip, shape);
163 
164  add(elvect, ip.weight*(Qvec*tangent), shape, elvect);
165  }
166 }
167 
169  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
170 {
171  int vdim = Q.GetVDim();
172  int dof = el.GetDof();
173 
174  double val,cf;
175 
176  shape.SetSize(dof); // vector of size dof
177 
178  elvect.SetSize(dof * vdim);
179  elvect = 0.0;
180 
181  const IntegrationRule *ir = IntRule;
182  if (ir == NULL)
183  {
184  int intorder = el.GetOrder() + 1;
185  ir = &IntRules.Get(el.GetGeomType(), intorder);
186  }
187 
188  for (int i = 0; i < ir->GetNPoints(); i++)
189  {
190  const IntegrationPoint &ip = ir->IntPoint(i);
191 
192  Tr.SetIntPoint (&ip);
193  val = Tr.Weight();
194 
195  el.CalcShape(ip, shape);
196  Q.Eval (Qvec, Tr, ip);
197 
198  for (int k = 0; k < vdim; k++)
199  {
200  cf = val * Qvec(k);
201 
202  for (int s = 0; s < dof; s++)
203  {
204  elvect(dof*k+s) += ip.weight * cf * shape(s);
205  }
206  }
207  }
208 }
209 
211  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
212 {
213  MFEM_ASSERT(vec_delta != NULL, "coefficient must be VectorDeltaCoefficient");
214  int vdim = Q.GetVDim();
215  int dof = fe.GetDof();
216 
217  shape.SetSize(dof);
218  fe.CalcPhysShape(Trans, shape);
219 
220  vec_delta->EvalDelta(Qvec, Trans, Trans.GetIntPoint());
221 
222  elvect.SetSize(dof*vdim);
223  DenseMatrix elvec_as_mat(elvect.GetData(), dof, vdim);
224  MultVWt(shape, Qvec, elvec_as_mat);
225 }
226 
227 
229  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
230 {
231  int vdim = Q.GetVDim();
232  int dof = el.GetDof();
233 
234  shape.SetSize(dof);
235  vec.SetSize(vdim);
236 
237  elvect.SetSize(dof * vdim);
238  elvect = 0.0;
239 
240  const IntegrationRule *ir = IntRule;
241  if (ir == NULL)
242  {
243  int intorder = el.GetOrder() + 1;
244  ir = &IntRules.Get(el.GetGeomType(), intorder);
245  }
246 
247  for (int i = 0; i < ir->GetNPoints(); i++)
248  {
249  const IntegrationPoint &ip = ir->IntPoint(i);
250 
251  Q.Eval(vec, Tr, ip);
252  Tr.SetIntPoint (&ip);
253  vec *= Tr.Weight() * ip.weight;
254  el.CalcShape(ip, shape);
255  for (int k = 0; k < vdim; k++)
256  for (int s = 0; s < dof; s++)
257  {
258  elvect(dof*k+s) += vec(k) * shape(s);
259  }
260  }
261 }
262 
264  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
265 {
266  int vdim = Q.GetVDim();
267  int dof = el.GetDof();
268 
269  shape.SetSize(dof);
270  vec.SetSize(vdim);
271 
272  elvect.SetSize(dof * vdim);
273  elvect = 0.0;
274 
275  const IntegrationRule *ir = IntRule;
276  if (ir == NULL)
277  {
278  int intorder = el.GetOrder() + 1;
279  ir = &IntRules.Get(Tr.FaceGeom, intorder);
280  }
281 
282  for (int i = 0; i < ir->GetNPoints(); i++)
283  {
284  const IntegrationPoint &ip = ir->IntPoint(i);
285  IntegrationPoint eip;
286  Tr.Loc1.Transform(ip, eip);
287 
288  Tr.Face->SetIntPoint(&ip);
289  Q.Eval(vec, *Tr.Face, ip);
290  vec *= Tr.Face->Weight() * ip.weight;
291  el.CalcShape(eip, shape);
292  for (int k = 0; k < vdim; k++)
293  {
294  for (int s = 0; s < dof; s++)
295  {
296  elvect(dof*k+s) += vec(k) * shape(s);
297  }
298  }
299  }
300 }
301 
302 
304  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
305 {
306  int dof = el.GetDof();
307  int spaceDim = Tr.GetSpaceDim();
308 
309  vshape.SetSize(dof,spaceDim);
310  vec.SetSize(spaceDim);
311 
312  elvect.SetSize(dof);
313  elvect = 0.0;
314 
315  const IntegrationRule *ir = IntRule;
316  if (ir == NULL)
317  {
318  // int intorder = 2*el.GetOrder() - 1; // ok for O(h^{k+1}) conv. in L2
319  int intorder = 2*el.GetOrder();
320  ir = &IntRules.Get(el.GetGeomType(), intorder);
321  }
322 
323  for (int i = 0; i < ir->GetNPoints(); i++)
324  {
325  const IntegrationPoint &ip = ir->IntPoint(i);
326 
327  Tr.SetIntPoint (&ip);
328  el.CalcVShape(Tr, vshape);
329 
330  QF.Eval (vec, Tr, ip);
331  vec *= ip.weight * Tr.Weight();
332 
333  vshape.AddMult (vec, elvect);
334  }
335 }
336 
338  const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
339 {
340  MFEM_ASSERT(vec_delta != NULL, "coefficient must be VectorDeltaCoefficient");
341  int dof = fe.GetDof();
342  int spaceDim = Trans.GetSpaceDim();
343 
344  vshape.SetSize(dof, spaceDim);
345  fe.CalcPhysVShape(Trans, vshape);
346 
347  vec_delta->EvalDelta(vec, Trans, Trans.GetIntPoint());
348 
349  elvect.SetSize(dof);
350  vshape.Mult(vec, elvect);
351 }
352 
353 
355  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
356 {
357  int dim = el.GetDim()+1;
358  int dof = el.GetDof();
359 
360  shape.SetSize (dof);
361  nor.SetSize (dim);
362  elvect.SetSize (dim*dof);
363 
364  const IntegrationRule *ir = IntRule;
365  if (ir == NULL)
366  {
367  ir = &IntRules.Get(el.GetGeomType(), el.GetOrder() + 1);
368  }
369 
370  elvect = 0.0;
371  for (int i = 0; i < ir->GetNPoints(); i++)
372  {
373  const IntegrationPoint &ip = ir->IntPoint(i);
374  Tr.SetIntPoint (&ip);
375  CalcOrtho(Tr.Jacobian(), nor);
376  el.CalcShape (ip, shape);
377  nor *= Sign * ip.weight * F -> Eval (Tr, ip);
378  for (int j = 0; j < dof; j++)
379  for (int k = 0; k < dim; k++)
380  {
381  elvect(dof*k+j) += nor(k) * shape(j);
382  }
383  }
384 }
385 
386 
388  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
389 {
390  int dof = el.GetDof();
391 
392  shape.SetSize(dof);
393  elvect.SetSize(dof);
394  elvect = 0.0;
395 
396  const IntegrationRule *ir = IntRule;
397  if (ir == NULL)
398  {
399  int intorder = 2*el.GetOrder(); // <----------
400  ir = &IntRules.Get(el.GetGeomType(), intorder);
401  }
402 
403  for (int i = 0; i < ir->GetNPoints(); i++)
404  {
405  const IntegrationPoint &ip = ir->IntPoint(i);
406 
407  Tr.SetIntPoint (&ip);
408  double val = ip.weight*F.Eval(Tr, ip);
409 
410  el.CalcShape(ip, shape);
411 
412  add(elvect, val, shape, elvect);
413  }
414 }
415 
416 
418  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
419 {
420  int dof = el.GetDof();
421  DenseMatrix vshape(dof, 2);
422  Vector f_loc(3);
423  Vector f_hat(2);
424 
425  elvect.SetSize(dof);
426  elvect = 0.0;
427 
428  const IntegrationRule *ir = IntRule;
429  if (ir == NULL)
430  {
431  int intorder = 2*el.GetOrder(); // <----------
432  ir = &IntRules.Get(el.GetGeomType(), intorder);
433  }
434 
435  for (int i = 0; i < ir->GetNPoints(); i++)
436  {
437  const IntegrationPoint &ip = ir->IntPoint(i);
438 
439  Tr.SetIntPoint(&ip);
440  f.Eval(f_loc, Tr, ip);
441  Tr.Jacobian().MultTranspose(f_loc, f_hat);
442  el.CalcVShape(ip, vshape);
443 
444  Swap<double>(f_hat(0), f_hat(1));
445  f_hat(0) = -f_hat(0);
446  f_hat *= ip.weight;
447  vshape.AddMult(f_hat, elvect);
448  }
449 }
450 
451 
453  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
454 {
455  mfem_error("BoundaryFlowIntegrator::AssembleRHSElementVect\n"
456  " is not implemented as boundary integrator!\n"
457  " Use LinearForm::AddBdrFaceIntegrator instead of\n"
458  " LinearForm::AddBoundaryIntegrator.");
459 }
460 
462  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
463 {
464  int dim, ndof, order;
465  double un, w, vu_data[3], nor_data[3];
466 
467  dim = el.GetDim();
468  ndof = el.GetDof();
469  Vector vu(vu_data, dim), nor(nor_data, dim);
470 
471  const IntegrationRule *ir = IntRule;
472  if (ir == NULL)
473  {
474  // Assuming order(u)==order(mesh)
475  order = Tr.Elem1->OrderW() + 2*el.GetOrder();
476  if (el.Space() == FunctionSpace::Pk)
477  {
478  order++;
479  }
480  ir = &IntRules.Get(Tr.FaceGeom, order);
481  }
482 
483  shape.SetSize(ndof);
484  elvect.SetSize(ndof);
485  elvect = 0.0;
486 
487  for (int p = 0; p < ir->GetNPoints(); p++)
488  {
489  const IntegrationPoint &ip = ir->IntPoint(p);
490  IntegrationPoint eip;
491  Tr.Loc1.Transform(ip, eip);
492  el.CalcShape(eip, shape);
493 
494  Tr.Face->SetIntPoint(&ip);
495 
496  u->Eval(vu, *Tr.Elem1, eip);
497 
498  if (dim == 1)
499  {
500  nor(0) = 2*eip.x - 1.0;
501  }
502  else
503  {
504  CalcOrtho(Tr.Face->Jacobian(), nor);
505  }
506 
507  un = vu * nor;
508  w = 0.5*alpha*un - beta*fabs(un);
509  w *= ip.weight*f->Eval(*Tr.Elem1, eip);
510  elvect.Add(w, shape);
511  }
512 }
513 
514 
516  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
517 {
518  mfem_error("DGDirichletLFIntegrator::AssembleRHSElementVect");
519 }
520 
522  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
523 {
524  int dim, ndof;
525  bool kappa_is_nonzero = (kappa != 0.);
526  double w;
527 
528  dim = el.GetDim();
529  ndof = el.GetDof();
530 
531  nor.SetSize(dim);
532  nh.SetSize(dim);
533  ni.SetSize(dim);
534  adjJ.SetSize(dim);
535  if (MQ)
536  {
537  mq.SetSize(dim);
538  }
539 
540  shape.SetSize(ndof);
541  dshape.SetSize(ndof, dim);
542  dshape_dn.SetSize(ndof);
543 
544  elvect.SetSize(ndof);
545  elvect = 0.0;
546 
547  const IntegrationRule *ir = IntRule;
548  if (ir == NULL)
549  {
550  // a simple choice for the integration order; is this OK?
551  int order = 2*el.GetOrder();
552  ir = &IntRules.Get(Tr.FaceGeom, order);
553  }
554 
555  for (int p = 0; p < ir->GetNPoints(); p++)
556  {
557  const IntegrationPoint &ip = ir->IntPoint(p);
558  IntegrationPoint eip;
559 
560  Tr.Loc1.Transform(ip, eip);
561  Tr.Face->SetIntPoint(&ip);
562  if (dim == 1)
563  {
564  nor(0) = 2*eip.x - 1.0;
565  }
566  else
567  {
568  CalcOrtho(Tr.Face->Jacobian(), nor);
569  }
570 
571  el.CalcShape(eip, shape);
572  el.CalcDShape(eip, dshape);
573  Tr.Elem1->SetIntPoint(&eip);
574  // compute uD through the face transformation
575  w = ip.weight * uD->Eval(*Tr.Face, ip) / Tr.Elem1->Weight();
576  if (!MQ)
577  {
578  if (Q)
579  {
580  w *= Q->Eval(*Tr.Elem1, eip);
581  }
582  ni.Set(w, nor);
583  }
584  else
585  {
586  nh.Set(w, nor);
587  MQ->Eval(mq, *Tr.Elem1, eip);
588  mq.MultTranspose(nh, ni);
589  }
591  adjJ.Mult(ni, nh);
592 
594  elvect.Add(sigma, dshape_dn);
595 
596  if (kappa_is_nonzero)
597  {
598  elvect.Add(kappa*(ni*nor), shape);
599  }
600  }
601 }
602 
603 
605  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
606 {
607  mfem_error("DGElasticityDirichletLFIntegrator::AssembleRHSElementVect");
608 }
609 
611  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
612 {
613  MFEM_ASSERT(Tr.Elem2No < 0, "interior boundary is not supported");
614 
615 #ifdef MFEM_THREAD_SAFE
616  Vector shape;
620  Vector nor;
623  Vector u_dir;
624 #endif
625 
626  const int dim = el.GetDim();
627  const int ndofs = el.GetDof();
628  const int nvdofs = dim*ndofs;
629 
630  elvect.SetSize(nvdofs);
631  elvect = 0.0;
632 
633  adjJ.SetSize(dim);
634  shape.SetSize(ndofs);
635  dshape.SetSize(ndofs, dim);
636  dshape_ps.SetSize(ndofs, dim);
637  nor.SetSize(dim);
638  dshape_dn.SetSize(ndofs);
639  dshape_du.SetSize(ndofs);
640  u_dir.SetSize(dim);
641 
642  const IntegrationRule *ir = IntRule;
643  if (ir == NULL)
644  {
645  const int order = 2*el.GetOrder(); // <-----
646  ir = &IntRules.Get(Tr.FaceGeom, order);
647  }
648 
649  for (int pi = 0; pi < ir->GetNPoints(); ++pi)
650  {
651  const IntegrationPoint &ip = ir->IntPoint(pi);
652  IntegrationPoint eip;
653  Tr.Loc1.Transform(ip, eip);
654  Tr.Face->SetIntPoint(&ip);
655  Tr.Elem1->SetIntPoint(&eip);
656 
657  // Evaluate the Dirichlet b.c. using the face transformation.
658  uD.Eval(u_dir, *Tr.Face, ip);
659 
660  el.CalcShape(eip, shape);
661  el.CalcDShape(eip, dshape);
662 
664  Mult(dshape, adjJ, dshape_ps);
665 
666  if (dim == 1)
667  {
668  nor(0) = 2*eip.x - 1.0;
669  }
670  else
671  {
672  CalcOrtho(Tr.Face->Jacobian(), nor);
673  }
674 
675  double wL, wM, jcoef;
676  {
677  const double w = ip.weight / Tr.Elem1->Weight();
678  wL = w * lambda->Eval(*Tr.Elem1, eip);
679  wM = w * mu->Eval(*Tr.Elem1, eip);
680  jcoef = kappa * (wL + 2.0*wM) * (nor*nor);
681  dshape_ps.Mult(nor, dshape_dn);
682  dshape_ps.Mult(u_dir, dshape_du);
683  }
684 
685  // alpha < uD, (lambda div(v) I + mu (grad(v) + grad(v)^T)) . n > +
686  // + kappa < h^{-1} (lambda + 2 mu) uD, v >
687 
688  // i = idof + ndofs * im
689  // v_phi(i,d) = delta(im,d) phi(idof)
690  // div(v_phi(i)) = dphi(idof,im)
691  // (grad(v_phi(i)))(k,l) = delta(im,k) dphi(idof,l)
692  //
693  // term 1:
694  // alpha < uD, lambda div(v_phi(i)) n >
695  // alpha lambda div(v_phi(i)) (uD.n) =
696  // alpha lambda dphi(idof,im) (uD.n) --> quadrature -->
697  // ip.weight/det(J1) alpha lambda (uD.nor) dshape_ps(idof,im) =
698  // alpha * wL * (u_dir*nor) * dshape_ps(idof,im)
699  // term 2:
700  // < alpha uD, mu grad(v_phi(i)).n > =
701  // alpha mu uD^T grad(v_phi(i)) n =
702  // alpha mu uD(k) delta(im,k) dphi(idof,l) n(l) =
703  // alpha mu uD(im) dphi(idof,l) n(l) --> quadrature -->
704  // ip.weight/det(J1) alpha mu uD(im) dshape_ps(idof,l) nor(l) =
705  // alpha * wM * u_dir(im) * dshape_dn(idof)
706  // term 3:
707  // < alpha uD, mu (grad(v_phi(i)))^T n > =
708  // alpha mu n^T grad(v_phi(i)) uD =
709  // alpha mu n(k) delta(im,k) dphi(idof,l) uD(l) =
710  // alpha mu n(im) dphi(idof,l) uD(l) --> quadrature -->
711  // ip.weight/det(J1) alpha mu nor(im) dshape_ps(idof,l) uD(l) =
712  // alpha * wM * nor(im) * dshape_du(idof)
713  // term j:
714  // < kappa h^{-1} (lambda + 2 mu) uD, v_phi(i) > =
715  // kappa/h (lambda + 2 mu) uD(k) v_phi(i,k) =
716  // kappa/h (lambda + 2 mu) uD(k) delta(im,k) phi(idof) =
717  // kappa/h (lambda + 2 mu) uD(im) phi(idof) --> quadrature -->
718  // [ 1/h = |nor|/det(J1) ]
719  // ip.weight/det(J1) |nor|^2 kappa (lambda + 2 mu) uD(im) phi(idof) =
720  // jcoef * u_dir(im) * shape(idof)
721 
722  wM *= alpha;
723  const double t1 = alpha * wL * (u_dir*nor);
724  for (int im = 0, i = 0; im < dim; ++im)
725  {
726  const double t2 = wM * u_dir(im);
727  const double t3 = wM * nor(im);
728  const double tj = jcoef * u_dir(im);
729  for (int idof = 0; idof < ndofs; ++idof, ++i)
730  {
731  elvect(i) += (t1*dshape_ps(idof,im) + t2*dshape_dn(idof) +
732  t3*dshape_du(idof) + tj*shape(idof));
733  }
734  }
735  }
736 }
737 
738 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:232
Abstract class for Finite Elements.
Definition: fe.hpp:140
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:387
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:211
ElementTransformation * Face
Definition: eltrans.hpp:347
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
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:3736
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:452
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:604
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:478
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
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)
Definition: eltrans.hpp:52
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:3042
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:221
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:168
int Space() const
Returns the type of space on each element.
Definition: fe.hpp:231
const IntegrationRule * IntRule
Definition: lininteg.hpp:25
const IntegrationPoint & GetIntPoint()
Definition: eltrans.hpp:54
Polynomials of order k.
Definition: fe.hpp:127
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:129
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:3272
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:67
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:260
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:235
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:219
int dim
Definition: ex3.cpp:47
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:348
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:146
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:67
int GetSpaceDim() const
Get the dimension of the target (physical) space.
Definition: eltrans.hpp:93
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:228
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...
VectorDeltaCoefficient * vec_delta
Definition: lininteg.hpp:51
MatrixCoefficient * MQ
Definition: lininteg.hpp:352
int GetVDim()
Returns dimension of the vector.
void AddMult(const Vector &x, Vector &y) const
y += A.x
Definition: densemat.cpp:242
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:128
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:515
void CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape) const
Equivalent to the CalcVShape() method with the same arguments.
Definition: fe.hpp:292
int GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:214
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:217
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:179
DeltaCoefficient * delta
Definition: lininteg.hpp:50
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:515
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:354
double EvalDelta(ElementTransformation &T, const IntegrationPoint &ip)
Return the Scale() multiplied by the weight Coefficient, if any.
Definition: coefficient.cpp:83
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:219
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:417
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:201
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:337
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:96
ElementTransformation * Elem1
Definition: eltrans.hpp:347
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Vector data type.
Definition: vector.hpp:48
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:168
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:303
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:210
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:86
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:353