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