MFEM  v3.2
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 &el, ElementTransformation &Tr, Vector &elvect)
59 {
60  int dof = el.GetDof();
61 
62  shape.SetSize(dof); // vector of size dof
63  elvect.SetSize(dof);
64  elvect = 0.0;
65 
66  const IntegrationRule *ir = IntRule;
67  if (ir == NULL)
68  {
69  int intorder = oa * el.GetOrder() + ob; // <----------
70  ir = &IntRules.Get(el.GetGeomType(), intorder);
71  }
72 
73  for (int i = 0; i < ir->GetNPoints(); i++)
74  {
75  const IntegrationPoint &ip = ir->IntPoint(i);
76 
77  Tr.SetIntPoint (&ip);
78  double val = Tr.Weight() * Q.Eval(Tr, ip);
79 
80  el.CalcShape(ip, shape);
81 
82  add(elvect, ip.weight * val, shape, elvect);
83  }
84 }
85 
87  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
88 {
89  int dim = el.GetDim()+1;
90  int dof = el.GetDof();
91  Vector nor(dim), Qvec;
92 
93  shape.SetSize(dof);
94  elvect.SetSize(dof);
95  elvect = 0.0;
96 
97  const IntegrationRule *ir = IntRule;
98  if (ir == NULL)
99  {
100  int intorder = oa * el.GetOrder() + ob; // <----------
101  ir = &IntRules.Get(el.GetGeomType(), intorder);
102  }
103 
104  for (int i = 0; i < ir->GetNPoints(); i++)
105  {
106  const IntegrationPoint &ip = ir->IntPoint(i);
107 
108  Tr.SetIntPoint(&ip);
109  CalcOrtho(Tr.Jacobian(), nor);
110  Q.Eval(Qvec, Tr, ip);
111 
112  el.CalcShape(ip, shape);
113 
114  elvect.Add(ip.weight*(Qvec*nor), shape);
115  }
116 }
117 
119  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
120 {
121  int dim = el.GetDim()+1;
122  int dof = el.GetDof();
123  Vector tangent(dim), Qvec;
124 
125  shape.SetSize(dof);
126  elvect.SetSize(dof);
127  elvect = 0.0;
128 
129  if (dim != 2)
130  {
131  mfem_error("These methods make sense only in 2D problems.");
132  }
133 
134  const IntegrationRule *ir = IntRule;
135  if (ir == NULL)
136  {
137  int intorder = oa * el.GetOrder() + ob; // <----------
138  ir = &IntRules.Get(el.GetGeomType(), intorder);
139  }
140 
141  for (int i = 0; i < ir->GetNPoints(); i++)
142  {
143  const IntegrationPoint &ip = ir->IntPoint(i);
144 
145  Tr.SetIntPoint(&ip);
146  const DenseMatrix &Jac = Tr.Jacobian();
147  tangent(0) = Jac(0,0);
148  tangent(1) = Jac(1,0);
149 
150  Q.Eval(Qvec, Tr, ip);
151 
152  el.CalcShape(ip, shape);
153 
154  add(elvect, ip.weight*(Qvec*tangent), shape, elvect);
155  }
156 }
157 
159  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
160 {
161  int vdim = Q.GetVDim();
162  int dof = el.GetDof();
163 
164  double val,cf;
165 
166  shape.SetSize(dof); // vector of size dof
167 
168  elvect.SetSize(dof * vdim);
169  elvect = 0.0;
170 
171  const IntegrationRule *ir = IntRule;
172  if (ir == NULL)
173  {
174  int intorder = el.GetOrder() + 1;
175  ir = &IntRules.Get(el.GetGeomType(), intorder);
176  }
177 
178  for (int i = 0; i < ir->GetNPoints(); i++)
179  {
180  const IntegrationPoint &ip = ir->IntPoint(i);
181 
182  Tr.SetIntPoint (&ip);
183  val = Tr.Weight();
184 
185  el.CalcShape(ip, shape);
186  Q.Eval (Qvec, Tr, ip);
187 
188  for (int k = 0; k < vdim; k++)
189  {
190  cf = val * Qvec(k);
191 
192  for (int s = 0; s < dof; s++)
193  {
194  elvect(dof*k+s) += ip.weight * cf * shape(s);
195  }
196  }
197  }
198 }
199 
201  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
202 {
203  int vdim = Q.GetVDim();
204  int dof = el.GetDof();
205 
206  shape.SetSize(dof);
207  vec.SetSize(vdim);
208 
209  elvect.SetSize(dof * vdim);
210  elvect = 0.0;
211 
212  const IntegrationRule *ir = IntRule;
213  if (ir == NULL)
214  {
215  int intorder = el.GetOrder() + 1;
216  ir = &IntRules.Get(el.GetGeomType(), intorder);
217  }
218 
219  for (int i = 0; i < ir->GetNPoints(); i++)
220  {
221  const IntegrationPoint &ip = ir->IntPoint(i);
222 
223  Q.Eval(vec, Tr, ip);
224  Tr.SetIntPoint (&ip);
225  vec *= Tr.Weight() * ip.weight;
226  el.CalcShape(ip, shape);
227  for (int k = 0; k < vdim; k++)
228  for (int s = 0; s < dof; s++)
229  {
230  elvect(dof*k+s) += vec(k) * shape(s);
231  }
232  }
233 }
234 
236  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
237 {
238  int dof = el.GetDof();
239  int spaceDim = Tr.GetSpaceDim();
240 
241  vshape.SetSize(dof,spaceDim);
242  vec.SetSize(spaceDim);
243 
244  elvect.SetSize(dof);
245  elvect = 0.0;
246 
247  const IntegrationRule *ir = IntRule;
248  if (ir == NULL)
249  {
250  // int intorder = 2*el.GetOrder() - 1; // ok for O(h^{k+1}) conv. in L2
251  int intorder = 2*el.GetOrder();
252  ir = &IntRules.Get(el.GetGeomType(), intorder);
253  }
254 
255  for (int i = 0; i < ir->GetNPoints(); i++)
256  {
257  const IntegrationPoint &ip = ir->IntPoint(i);
258 
259  Tr.SetIntPoint (&ip);
260  el.CalcVShape(Tr, vshape);
261 
262  QF.Eval (vec, Tr, ip);
263  vec *= ip.weight * Tr.Weight();
264 
265  vshape.AddMult (vec, elvect);
266  }
267 }
268 
269 
271  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
272 {
273  int dim = el.GetDim()+1;
274  int dof = el.GetDof();
275 
276  shape.SetSize (dof);
277  nor.SetSize (dim);
278  elvect.SetSize (dim*dof);
279 
280  const IntegrationRule *ir = IntRule;
281  if (ir == NULL)
282  {
283  ir = &IntRules.Get(el.GetGeomType(), el.GetOrder() + 1);
284  }
285 
286  elvect = 0.0;
287  for (int i = 0; i < ir->GetNPoints(); i++)
288  {
289  const IntegrationPoint &ip = ir->IntPoint(i);
290  Tr.SetIntPoint (&ip);
291  CalcOrtho(Tr.Jacobian(), nor);
292  el.CalcShape (ip, shape);
293  nor *= Sign * ip.weight * F -> Eval (Tr, ip);
294  for (int j = 0; j < dof; j++)
295  for (int k = 0; k < dim; k++)
296  {
297  elvect(dof*k+j) += nor(k) * shape(j);
298  }
299  }
300 }
301 
302 
304  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
305 {
306  int dof = el.GetDof();
307 
308  shape.SetSize(dof);
309  elvect.SetSize(dof);
310  elvect = 0.0;
311 
312  const IntegrationRule *ir = IntRule;
313  if (ir == NULL)
314  {
315  int intorder = 2*el.GetOrder(); // <----------
316  ir = &IntRules.Get(el.GetGeomType(), intorder);
317  }
318 
319  for (int i = 0; i < ir->GetNPoints(); i++)
320  {
321  const IntegrationPoint &ip = ir->IntPoint(i);
322 
323  Tr.SetIntPoint (&ip);
324  double val = ip.weight*F.Eval(Tr, ip);
325 
326  el.CalcShape(ip, shape);
327 
328  add(elvect, val, shape, elvect);
329  }
330 }
331 
332 
334  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
335 {
336  int dof = el.GetDof();
337  DenseMatrix vshape(dof, 2);
338  Vector f_loc(3);
339  Vector f_hat(2);
340 
341  elvect.SetSize(dof);
342  elvect = 0.0;
343 
344  const IntegrationRule *ir = IntRule;
345  if (ir == NULL)
346  {
347  int intorder = 2*el.GetOrder(); // <----------
348  ir = &IntRules.Get(el.GetGeomType(), intorder);
349  }
350 
351  for (int i = 0; i < ir->GetNPoints(); i++)
352  {
353  const IntegrationPoint &ip = ir->IntPoint(i);
354 
355  Tr.SetIntPoint(&ip);
356  f.Eval(f_loc, Tr, ip);
357  Tr.Jacobian().MultTranspose(f_loc, f_hat);
358  el.CalcVShape(ip, vshape);
359 
360  Swap<double>(f_hat(0), f_hat(1));
361  f_hat(0) = -f_hat(0);
362  f_hat *= ip.weight;
363  vshape.AddMult(f_hat, elvect);
364  }
365 }
366 
367 
369  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
370 {
371  mfem_error("BoundaryFlowIntegrator::AssembleRHSElementVect\n"
372  " is not implemented as boundary integrator!\n"
373  " Use LinearForm::AddBdrFaceIntegrator instead of\n"
374  " LinearForm::AddBoundaryIntegrator.");
375 }
376 
378  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
379 {
380  int dim, ndof, order;
381  double un, w, vu_data[3], nor_data[3];
382 
383  dim = el.GetDim();
384  ndof = el.GetDof();
385  Vector vu(vu_data, dim), nor(nor_data, dim);
386 
387  const IntegrationRule *ir = IntRule;
388  if (ir == NULL)
389  {
390  // Assuming order(u)==order(mesh)
391  order = Tr.Elem1->OrderW() + 2*el.GetOrder();
392  if (el.Space() == FunctionSpace::Pk)
393  {
394  order++;
395  }
396  ir = &IntRules.Get(Tr.FaceGeom, order);
397  }
398 
399  shape.SetSize(ndof);
400  elvect.SetSize(ndof);
401  elvect = 0.0;
402 
403  for (int p = 0; p < ir->GetNPoints(); p++)
404  {
405  const IntegrationPoint &ip = ir->IntPoint(p);
406  IntegrationPoint eip;
407  Tr.Loc1.Transform(ip, eip);
408  el.CalcShape(eip, shape);
409 
410  Tr.Face->SetIntPoint(&ip);
411 
412  u->Eval(vu, *Tr.Elem1, eip);
413 
414  if (dim == 1)
415  {
416  nor(0) = 2*eip.x - 1.0;
417  }
418  else
419  {
420  CalcOrtho(Tr.Face->Jacobian(), nor);
421  }
422 
423  un = vu * nor;
424  w = 0.5*alpha*un - beta*fabs(un);
425  w *= ip.weight*f->Eval(*Tr.Elem1, eip);
426  elvect.Add(w, shape);
427  }
428 }
429 
430 
432  const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
433 {
434  mfem_error("DGDirichletLFIntegrator::AssembleRHSElementVect");
435 }
436 
438  const FiniteElement &el, FaceElementTransformations &Tr, Vector &elvect)
439 {
440  int dim, ndof;
441  bool kappa_is_nonzero = (kappa != 0.);
442  double w;
443 
444  dim = el.GetDim();
445  ndof = el.GetDof();
446 
447  nor.SetSize(dim);
448  nh.SetSize(dim);
449  ni.SetSize(dim);
450  adjJ.SetSize(dim);
451  if (MQ)
452  {
453  mq.SetSize(dim);
454  }
455 
456  shape.SetSize(ndof);
457  dshape.SetSize(ndof, dim);
458  dshape_dn.SetSize(ndof);
459 
460  elvect.SetSize(ndof);
461  elvect = 0.0;
462 
463  const IntegrationRule *ir = IntRule;
464  if (ir == NULL)
465  {
466  // a simple choice for the integration order; is this OK?
467  int order = 2*el.GetOrder();
468  ir = &IntRules.Get(Tr.FaceGeom, order);
469  }
470 
471  for (int p = 0; p < ir->GetNPoints(); p++)
472  {
473  const IntegrationPoint &ip = ir->IntPoint(p);
474  IntegrationPoint eip;
475 
476  Tr.Loc1.Transform(ip, eip);
477  Tr.Face->SetIntPoint(&ip);
478  if (dim == 1)
479  {
480  nor(0) = 2*eip.x - 1.0;
481  }
482  else
483  {
484  CalcOrtho(Tr.Face->Jacobian(), nor);
485  }
486 
487  el.CalcShape(eip, shape);
488  el.CalcDShape(eip, dshape);
489  Tr.Elem1->SetIntPoint(&eip);
490  // compute uD through the face transformation
491  w = ip.weight * uD->Eval(*Tr.Face, ip) / Tr.Elem1->Weight();
492  if (!MQ)
493  {
494  if (Q)
495  {
496  w *= Q->Eval(*Tr.Elem1, eip);
497  }
498  ni.Set(w, nor);
499  }
500  else
501  {
502  nh.Set(w, nor);
503  MQ->Eval(mq, *Tr.Elem1, eip);
504  mq.MultTranspose(nh, ni);
505  }
507  adjJ.Mult(ni, nh);
508 
510  elvect.Add(sigma, dshape_dn);
511 
512  if (kappa_is_nonzero)
513  {
514  elvect.Add(kappa*(ni*nor), shape);
515  }
516  }
517 }
518 
519 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:228
Abstract class for Finite Elements.
Definition: fe.hpp:44
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:303
int GetDim() const
Returns the space dimension for the finite element.
Definition: fe.hpp:82
ElementTransformation * Face
Definition: eltrans.hpp:123
Class for integration rule.
Definition: intrules.hpp:83
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Definition: fe.cpp:30
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:235
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
void SetSize(int s)
Resize the vector if the new size is different.
Definition: vector.hpp:263
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:368
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:35
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2904
int GetOrder() const
Returns the order of the finite element.
Definition: fe.hpp:91
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:158
int Space() const
Returns the type of space on each element.
Definition: fe.hpp:94
const IntegrationRule * IntRule
Definition: lininteg.hpp:25
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:3134
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:57
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:259
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:231
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:193
int dim
Definition: ex3.cpp:47
IntegrationRules IntRules(0)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:292
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:124
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:200
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
MatrixCoefficient * MQ
Definition: lininteg.hpp:287
virtual double Weight()=0
int GetVDim()
Returns dimension of the vector.
void AddMult(const Vector &x, Vector &y) const
y += A.x
Definition: densemat.cpp:216
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:118
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:431
int GetGeomType() const
Returns the geometry type:
Definition: fe.hpp:85
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)=0
int GetDof() const
Returns the degrees of freedom in the FE space.
Definition: fe.hpp:88
void mfem_error(const char *msg)
Definition: error.cpp:23
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:243
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:270
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:218
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:333
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:200
Class for integration point with weight.
Definition: intrules.hpp:25
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:86
ElementTransformation * Elem1
Definition: eltrans.hpp:123
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Vector data type.
Definition: vector.hpp:33
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:142
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: lininteg.cpp:235
virtual int GetSpaceDim()=0
virtual const DenseMatrix & Jacobian()=0
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:71
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0