MFEM  v3.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
bilininteg.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 // Implementation of Bilinear Form Integrators
13 
14 #include "fem.hpp"
15 #include <cmath>
16 #include <algorithm>
17 
18 using namespace std;
19 
20 namespace mfem
21 {
22 
23 void BilinearFormIntegrator::AssembleElementMatrix (
24  const FiniteElement &el, ElementTransformation &Trans,
25  DenseMatrix &elmat )
26 {
27  mfem_error ("BilinearFormIntegrator::AssembleElementMatrix (...)\n"
28  " is not implemented fot this class.");
29 }
30 
31 void BilinearFormIntegrator::AssembleElementMatrix2 (
32  const FiniteElement &el1, const FiniteElement &el2,
33  ElementTransformation &Trans, DenseMatrix &elmat )
34 {
35  mfem_error ("BilinearFormIntegrator::AssembleElementMatrix2 (...)\n"
36  " is not implemented fot this class.");
37 }
38 
39 void BilinearFormIntegrator::AssembleFaceMatrix (
40  const FiniteElement &el1, const FiniteElement &el2,
42 {
43  mfem_error ("BilinearFormIntegrator::AssembleFaceMatrix (...)\n"
44  " is not implemented fot this class.");
45 }
46 
47 void BilinearFormIntegrator::AssembleFaceMatrix(
48  const FiniteElement &trial_face_fe, const FiniteElement &test_fe1,
49  const FiniteElement &test_fe2, FaceElementTransformations &Trans,
50  DenseMatrix &elmat)
51 {
52  MFEM_ABORT("AssembleFaceMatrix (mixed form) is not implemented for this"
53  " Integrator class.");
54 }
55 
56 void BilinearFormIntegrator::AssembleElementVector(
57  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
58  Vector &elvect)
59 {
60  mfem_error("BilinearFormIntegrator::AssembleElementVector\n"
61  " is not implemented fot this class.");
62 }
63 
64 
65 void TransposeIntegrator::AssembleElementMatrix (
66  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
67 {
68  bfi -> AssembleElementMatrix (el, Trans, bfi_elmat);
69  // elmat = bfi_elmat^t
70  elmat.Transpose (bfi_elmat);
71 }
72 
73 void TransposeIntegrator::AssembleElementMatrix2 (
74  const FiniteElement &trial_fe, const FiniteElement &test_fe,
75  ElementTransformation &Trans, DenseMatrix &elmat)
76 {
77  bfi -> AssembleElementMatrix2 (test_fe, trial_fe, Trans, bfi_elmat);
78  // elmat = bfi_elmat^t
79  elmat.Transpose (bfi_elmat);
80 }
81 
82 void TransposeIntegrator::AssembleFaceMatrix (
83  const FiniteElement &el1, const FiniteElement &el2,
85 {
86  bfi -> AssembleFaceMatrix (el1, el2, Trans, bfi_elmat);
87  // elmat = bfi_elmat^t
88  elmat.Transpose (bfi_elmat);
89 }
90 
91 void LumpedIntegrator::AssembleElementMatrix (
92  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
93 {
94  bfi -> AssembleElementMatrix (el, Trans, elmat);
95  elmat.Lump();
96 }
97 
98 void InverseIntegrator::AssembleElementMatrix(
99  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
100 {
101  integrator->AssembleElementMatrix(el, Trans, elmat);
102  elmat.Invert();
103 }
104 
105 void SumIntegrator::AssembleElementMatrix(
106  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
107 {
108  MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
109 
110  integrators[0]->AssembleElementMatrix(el, Trans, elmat);
111  for (int i = 1; i < integrators.Size(); i++)
112  {
113  integrators[i]->AssembleElementMatrix(el, Trans, elem_mat);
114  elmat += elem_mat;
115  }
116 }
117 
118 SumIntegrator::~SumIntegrator()
119 {
120  if (own_integrators)
121  {
122  for (int i = 0; i < integrators.Size(); i++)
123  delete integrators[i];
124  }
125 }
126 
127 
128 void DiffusionIntegrator::AssembleElementMatrix
130  DenseMatrix &elmat )
131 {
132  int nd = el.GetDof();
133  int dim = el.GetDim();
134  int spaceDim = Trans.GetSpaceDim();
135  bool square = (dim == spaceDim);
136  double w;
137 
138 #ifdef MFEM_THREAD_SAFE
139  DenseMatrix dshape(nd,dim), dshapedxt(nd,spaceDim), invdfdx(dim,spaceDim);
140 #else
141  dshape.SetSize(nd,dim);
142  dshapedxt.SetSize(nd,spaceDim);
143  invdfdx.SetSize(dim,spaceDim);
144 #endif
145  elmat.SetSize(nd);
146 
147  const IntegrationRule *ir = IntRule;
148  if (ir == NULL)
149  {
150  int order;
151  if (el.Space() == FunctionSpace::Pk)
152  order = 2*el.GetOrder() - 2;
153  else
154  // order = 2*el.GetOrder() - 2; // <-- this seems to work fine too
155  order = 2*el.GetOrder() + dim - 1;
156 
157  if (el.Space() == FunctionSpace::rQk)
158  ir = &RefinedIntRules.Get(el.GetGeomType(), order);
159  else
160  ir = &IntRules.Get(el.GetGeomType(), order);
161  }
162 
163  elmat = 0.0;
164  for (int i = 0; i < ir->GetNPoints(); i++)
165  {
166  const IntegrationPoint &ip = ir->IntPoint(i);
167  el.CalcDShape(ip, dshape);
168 
169  Trans.SetIntPoint(&ip);
170  // Compute invdfdx = / adj(J), if J is square
171  // \ adj(J^t.J).J^t, otherwise
172  CalcAdjugate(Trans.Jacobian(), invdfdx);
173  w = Trans.Weight();
174  w = ip.weight / (square ? w : w*w*w);
175  Mult(dshape, invdfdx, dshapedxt);
176  if (!MQ)
177  {
178  if (Q)
179  w *= Q->Eval(Trans, ip);
180  AddMult_a_AAt(w, dshapedxt, elmat);
181  }
182  else
183  {
184  MQ->Eval(invdfdx, Trans, ip);
185  invdfdx *= w;
186  Mult(dshapedxt, invdfdx, dshape);
187  AddMultABt(dshape, dshapedxt, elmat);
188  }
189  }
190 }
191 
192 void DiffusionIntegrator::AssembleElementMatrix2(
193  const FiniteElement &trial_fe, const FiniteElement &test_fe,
194  ElementTransformation &Trans, DenseMatrix &elmat)
195 {
196  int tr_nd = trial_fe.GetDof();
197  int te_nd = test_fe.GetDof();
198  int dim = trial_fe.GetDim();
199  int spaceDim = Trans.GetSpaceDim();
200  bool square = (dim == spaceDim);
201  double w;
202 
203 #ifdef MFEM_THREAD_SAFE
204  DenseMatrix dshape(tr_nd, dim), dshapedxt(tr_nd, spaceDim);
205  DenseMatrix te_dshape(te_nd, dim), te_dshapedxt(te_nd, spaceDim);
206  DenseMatrix invdfdx(dim, spaceDim);
207 #else
208  dshape.SetSize(tr_nd, dim);
209  dshapedxt.SetSize(tr_nd, spaceDim);
210  te_dshape.SetSize(te_nd, dim);
211  te_dshapedxt.SetSize(te_nd, spaceDim);
212  invdfdx.SetSize(dim, spaceDim);
213 #endif
214  elmat.SetSize(te_nd, tr_nd);
215 
216  const IntegrationRule *ir = IntRule;
217  if (ir == NULL)
218  {
219  int order;
220  if (trial_fe.Space() == FunctionSpace::Pk)
221  order = trial_fe.GetOrder() + test_fe.GetOrder() - 2;
222  else
223  order = trial_fe.GetOrder() + test_fe.GetOrder() + dim - 1;
224 
225  if (trial_fe.Space() == FunctionSpace::rQk)
226  ir = &RefinedIntRules.Get(trial_fe.GetGeomType(), order);
227  else
228  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
229  }
230 
231  elmat = 0.0;
232  for (int i = 0; i < ir->GetNPoints(); i++)
233  {
234  const IntegrationPoint &ip = ir->IntPoint(i);
235  trial_fe.CalcDShape(ip, dshape);
236  test_fe.CalcDShape(ip, te_dshape);
237 
238  Trans.SetIntPoint(&ip);
239  CalcAdjugate(Trans.Jacobian(), invdfdx);
240  w = Trans.Weight();
241  w = ip.weight / (square ? w : w*w*w);
242  Mult(dshape, invdfdx, dshapedxt);
243  Mult(te_dshape, invdfdx, te_dshapedxt);
244  // invdfdx, dshape, and te_dshape no longer needed
245  if (!MQ)
246  {
247  if (Q)
248  w *= Q->Eval(Trans, ip);
249  dshapedxt *= w;
250  AddMultABt(te_dshapedxt, dshapedxt, elmat);
251  }
252  else
253  {
254  MQ->Eval(invdfdx, Trans, ip);
255  invdfdx *= w;
256  Mult(te_dshapedxt, invdfdx, te_dshape);
257  AddMultABt(te_dshape, dshapedxt, elmat);
258  }
259  }
260 }
261 
262 void DiffusionIntegrator::AssembleElementVector(
263  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
264  Vector &elvect)
265 {
266  int nd = el.GetDof();
267  int dim = el.GetDim();
268  double w;
269 
270 #ifdef MFEM_THREAD_SAFE
271  DenseMatrix dshape(nd,dim), invdfdx(dim), mq(dim);
272 #else
273  dshape.SetSize(nd,dim);
274  invdfdx.SetSize(dim);
275  mq.SetSize(dim);
276 #endif
277  vec.SetSize(dim);
278  pointflux.SetSize(dim);
279 
280  elvect.SetSize(nd);
281 
282  const IntegrationRule *ir = IntRule;
283  if (ir == NULL)
284  {
285  int order;
286  if (el.Space() == FunctionSpace::Pk)
287  order = 2*el.GetOrder() - 2;
288  else
289  // order = 2*el.GetOrder() - 2; // <-- this seems to work fine too
290  order = 2*el.GetOrder() + dim - 1;
291 
292  if (el.Space() == FunctionSpace::rQk)
293  ir = &RefinedIntRules.Get(el.GetGeomType(), order);
294  else
295  ir = &IntRules.Get(el.GetGeomType(), order);
296  }
297 
298  elvect = 0.0;
299  for (int i = 0; i < ir->GetNPoints(); i++)
300  {
301  const IntegrationPoint &ip = ir->IntPoint(i);
302  el.CalcDShape(ip, dshape);
303 
304  Tr.SetIntPoint(&ip);
305  CalcAdjugate(Tr.Jacobian(), invdfdx); // invdfdx = adj(J)
306  w = ip.weight / Tr.Weight();
307 
308  if (!MQ)
309  {
310  dshape.MultTranspose(elfun, vec);
311  invdfdx.MultTranspose(vec, pointflux);
312  if (Q)
313  w *= Q->Eval(Tr, ip);
314  }
315  else
316  {
317 
318  dshape.MultTranspose(elfun, pointflux);
319  invdfdx.MultTranspose(pointflux, vec);
320  MQ->Eval(mq, Tr, ip);
321  mq.Mult(vec, pointflux);
322  }
323  pointflux *= w;
324  invdfdx.Mult(pointflux, vec);
325  dshape.AddMult(vec, elvect);
326  }
327 }
328 
329 void DiffusionIntegrator::ComputeElementFlux
331  Vector &u, const FiniteElement &fluxelem, Vector &flux, int wcoef )
332 {
333  int i, j, nd, dim, spaceDim, fnd;
334 
335  nd = el.GetDof();
336  dim = el.GetDim();
337  spaceDim = Trans.GetSpaceDim();
338 
339 #ifdef MFEM_THREAD_SAFE
340  DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim);
341 #else
342  dshape.SetSize(nd,dim);
343  invdfdx.SetSize(dim, spaceDim);
344 #endif
345  vec.SetSize(dim);
346  pointflux.SetSize(spaceDim);
347 
348  const IntegrationRule &ir = fluxelem.GetNodes();
349  fnd = ir.GetNPoints();
350  flux.SetSize( fnd * dim );
351 
352  for (i = 0; i < fnd; i++)
353  {
354  const IntegrationPoint &ip = ir.IntPoint(i);
355  el.CalcDShape(ip, dshape);
356  dshape.MultTranspose(u, vec);
357 
358  Trans.SetIntPoint (&ip);
359  CalcInverse(Trans.Jacobian(), invdfdx);
360  invdfdx.MultTranspose(vec, pointflux);
361 
362  if (!wcoef)
363  {
364  for (j = 0; j < dim; j++)
365  flux(fnd*j+i) = pointflux(j);
366  }
367  else if (!MQ)
368  {
369  if (Q)
370  pointflux *= Q->Eval(Trans,ip);
371  for (j = 0; j < dim; j++)
372  flux(fnd*j+i) = pointflux(j);
373  }
374  else
375  {
376  MQ->Eval(invdfdx, Trans, ip);
377  invdfdx.Mult(pointflux, vec);
378  for (j = 0; j < dim; j++)
379  flux(fnd*j+i) = vec(j);
380  }
381  }
382 }
383 
384 double DiffusionIntegrator::ComputeFluxEnergy
385 ( const FiniteElement &fluxelem, ElementTransformation &Trans,
386  Vector &flux)
387 {
388  int i, j, k, nd, dim, order;
389  double energy, co;
390 
391  nd = fluxelem.GetDof();
392  dim = fluxelem.GetDim();
393 
394 #ifdef MFEM_THREAD_SAFE
395  DenseMatrix invdfdx;
396 #endif
397 
398  shape.SetSize(nd);
399  pointflux.SetSize(dim);
400  if (MQ)
401  {
402  invdfdx.SetSize(dim);
403  vec.SetSize(dim);
404  }
405 
406  order = 2 * fluxelem.GetOrder(); // <--
407  const IntegrationRule *ir = &IntRules.Get(fluxelem.GetGeomType(), order);
408 
409  energy = 0.0;
410  for (i = 0; i < ir->GetNPoints(); i++)
411  {
412  const IntegrationPoint &ip = ir->IntPoint(i);
413  fluxelem.CalcShape(ip, shape);
414 
415  pointflux = 0.0;
416  for (k = 0; k < dim; k++)
417  for (j = 0; j < nd; j++)
418  pointflux(k) += flux(k*nd+j)*shape(j);
419 
420  Trans.SetIntPoint (&ip);
421  co = Trans.Weight() * ip.weight;
422 
423  if (!MQ)
424  {
425  co *= ( pointflux * pointflux );
426  if (Q)
427  co *= Q->Eval(Trans, ip);
428  }
429  else
430  {
431  MQ->Eval(invdfdx, Trans, ip);
432  co *= invdfdx.InnerProduct(pointflux, pointflux);
433  }
434 
435  energy += co;
436  }
437 
438  return energy;
439 }
440 
441 
442 void MassIntegrator::AssembleElementMatrix
444  DenseMatrix &elmat )
445 {
446  int nd = el.GetDof();
447  // int dim = el.GetDim();
448  double w;
449 
450  elmat.SetSize(nd);
451  shape.SetSize(nd);
452 
453  const IntegrationRule *ir = IntRule;
454  if (ir == NULL)
455  {
456  // int order = 2 * el.GetOrder();
457  int order = 2 * el.GetOrder() + Trans.OrderW();
458 
459  if (el.Space() == FunctionSpace::rQk)
460  ir = &RefinedIntRules.Get(el.GetGeomType(), order);
461  else
462  ir = &IntRules.Get(el.GetGeomType(), order);
463  }
464 
465  elmat = 0.0;
466  for (int i = 0; i < ir->GetNPoints(); i++)
467  {
468  const IntegrationPoint &ip = ir->IntPoint(i);
469  el.CalcShape(ip, shape);
470 
471  Trans.SetIntPoint (&ip);
472  w = Trans.Weight() * ip.weight;
473  if (Q)
474  w *= Q -> Eval(Trans, ip);
475 
476  AddMult_a_VVt(w, shape, elmat);
477  }
478 }
479 
480 void MassIntegrator::AssembleElementMatrix2(
481  const FiniteElement &trial_fe, const FiniteElement &test_fe,
482  ElementTransformation &Trans, DenseMatrix &elmat)
483 {
484  int tr_nd = trial_fe.GetDof();
485  int te_nd = test_fe.GetDof();
486  // int dim = trial_fe.GetDim();
487  double w;
488 
489  elmat.SetSize (te_nd, tr_nd);
490  shape.SetSize (tr_nd);
491  te_shape.SetSize (te_nd);
492 
493  const IntegrationRule *ir = IntRule;
494  if (ir == NULL)
495  {
496  int order = trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW();
497 
498  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
499  }
500 
501  elmat = 0.0;
502  for (int i = 0; i < ir->GetNPoints(); i++)
503  {
504  const IntegrationPoint &ip = ir->IntPoint(i);
505  trial_fe.CalcShape(ip, shape);
506  test_fe.CalcShape(ip, te_shape);
507 
508  Trans.SetIntPoint (&ip);
509  w = Trans.Weight() * ip.weight;
510  if (Q)
511  w *= Q -> Eval(Trans, ip);
512 
513  te_shape *= w;
514  AddMultVWt(te_shape, shape, elmat);
515  }
516 }
517 
518 
519 void ConvectionIntegrator::AssembleElementMatrix(
520  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
521 {
522  int nd = el.GetDof();
523  int dim = el.GetDim();
524 
525  elmat.SetSize(nd);
526  dshape.SetSize(nd,dim);
527  adjJ.SetSize(dim);
528  shape.SetSize(nd);
529  vec2.SetSize(dim);
530  BdFidxT.SetSize(nd);
531 
532  Vector vec1;
533 
534  const IntegrationRule *ir = IntRule;
535  if (ir == NULL)
536  {
537  int order = Trans.OrderGrad(&el) + Trans.Order() + el.GetOrder();
538  ir = &IntRules.Get(el.GetGeomType(), order);
539  }
540 
541  Q.Eval(Q_ir, Trans, *ir);
542 
543  elmat = 0.0;
544  for (int i = 0; i < ir->GetNPoints(); i++)
545  {
546  const IntegrationPoint &ip = ir->IntPoint(i);
547  el.CalcDShape(ip, dshape);
548  el.CalcShape(ip, shape);
549 
550  Trans.SetIntPoint(&ip);
551  CalcAdjugate(Trans.Jacobian(), adjJ);
552  Q_ir.GetColumnReference(i, vec1);
553  vec1 *= alpha * ip.weight;
554 
555  adjJ.Mult(vec1, vec2);
556  dshape.Mult(vec2, BdFidxT);
557 
558  AddMultVWt(shape, BdFidxT, elmat);
559  }
560 }
561 
562 
563 void GroupConvectionIntegrator::AssembleElementMatrix(
564  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
565 {
566  int nd = el.GetDof();
567  int dim = el.GetDim();
568 
569  elmat.SetSize(nd);
570  dshape.SetSize(nd,dim);
571  adjJ.SetSize(dim);
572  shape.SetSize(nd);
573  grad.SetSize(nd,dim);
574 
575  const IntegrationRule *ir = IntRule;
576  if (ir == NULL)
577  {
578  int order = Trans.OrderGrad(&el) + el.GetOrder();
579  ir = &IntRules.Get(el.GetGeomType(), order);
580  }
581 
582  Q.Eval(Q_nodal, Trans, el.GetNodes()); // sets the size of Q_nodal
583 
584  elmat = 0.0;
585  for (int i = 0; i < ir->GetNPoints(); i++)
586  {
587  const IntegrationPoint &ip = ir->IntPoint(i);
588  el.CalcDShape(ip, dshape);
589  el.CalcShape(ip, shape);
590 
591  Trans.SetIntPoint(&ip);
592  CalcAdjugate(Trans.Jacobian(), adjJ);
593 
594  Mult(dshape, adjJ, grad);
595 
596  double w = alpha * ip.weight;
597 
598  // elmat(k,l) += \sum_s w*shape(k)*Q_nodal(s,k)*grad(l,s)
599  for (int k = 0; k < nd; k++)
600  {
601  double wsk = w*shape(k);
602  for (int l = 0; l < nd; l++)
603  {
604  double a = 0.0;
605  for (int s = 0; s < dim; s++)
606  a += Q_nodal(s,k)*grad(l,s);
607  elmat(k,l) += wsk*a;
608  }
609  }
610  }
611 }
612 
613 
614 void VectorMassIntegrator::AssembleElementMatrix
616  DenseMatrix &elmat )
617 {
618  int nd = el.GetDof();
619  int dim = el.GetDim();
620  int vdim;
621 
622  double norm;
623 
624  // Get vdim from the ElementTransformation Trans ?
625  vdim = (VQ) ? (VQ -> GetVDim()) : ((MQ) ? (MQ -> GetVDim()) : (dim));
626 
627  elmat.SetSize(nd*vdim);
628  shape.SetSize(nd);
629  partelmat.SetSize(nd);
630  if (VQ)
631  vec.SetSize(vdim);
632  else if (MQ)
633  mcoeff.SetSize(vdim);
634 
635  const IntegrationRule *ir = IntRule;
636  if (ir == NULL)
637  {
638  int order = 2 * el.GetOrder() + Trans.OrderW() + Q_order;
639 
640  if (el.Space() == FunctionSpace::rQk)
641  ir = &RefinedIntRules.Get(el.GetGeomType(), order);
642  else
643  ir = &IntRules.Get(el.GetGeomType(), order);
644  }
645 
646  elmat = 0.0;
647  for (int s = 0; s < ir->GetNPoints(); s++)
648  {
649  const IntegrationPoint &ip = ir->IntPoint(s);
650  el.CalcShape(ip, shape);
651 
652  Trans.SetIntPoint (&ip);
653  norm = ip.weight * Trans.Weight();
654 
655  MultVVt(shape, partelmat);
656 
657  if (VQ)
658  {
659  VQ->Eval(vec, Trans, ip);
660  for (int k = 0; k < vdim; k++)
661  elmat.AddMatrix(norm*vec(k), partelmat, nd*k, nd*k);
662  }
663  else if (MQ)
664  {
665  MQ->Eval(mcoeff, Trans, ip);
666  for (int i = 0; i < vdim; i++)
667  for (int j = 0; j < vdim; j++)
668  elmat.AddMatrix(norm*mcoeff(i,j), partelmat, nd*i, nd*j);
669  }
670  else
671  {
672  if (Q)
673  norm *= Q->Eval(Trans, ip);
674  partelmat *= norm;
675  for (int k = 0; k < vdim; k++)
676  elmat.AddMatrix(partelmat, nd*k, nd*k);
677  }
678  }
679 }
680 
681 void VectorMassIntegrator::AssembleElementMatrix2(
682  const FiniteElement &trial_fe, const FiniteElement &test_fe,
683  ElementTransformation &Trans, DenseMatrix &elmat)
684 {
685  int tr_nd = trial_fe.GetDof();
686  int te_nd = test_fe.GetDof();
687  int dim = trial_fe.GetDim();
688  int vdim;
689 
690  double norm;
691 
692  // Get vdim from the ElementTransformation Trans ?
693  vdim = (VQ) ? (VQ -> GetVDim()) : ((MQ) ? (MQ -> GetVDim()) : (dim));
694 
695  elmat.SetSize(te_nd*vdim, tr_nd*vdim);
696  shape.SetSize(tr_nd);
697  te_shape.SetSize(te_nd);
698  partelmat.SetSize(te_nd, tr_nd);
699  if (VQ)
700  vec.SetSize(vdim);
701  else if (MQ)
702  mcoeff.SetSize(vdim);
703 
704  const IntegrationRule *ir = IntRule;
705  if (ir == NULL)
706  {
707  int order = (trial_fe.GetOrder() + test_fe.GetOrder() +
708  Trans.OrderW() + Q_order);
709 
710  if (trial_fe.Space() == FunctionSpace::rQk)
711  ir = &RefinedIntRules.Get(trial_fe.GetGeomType(), order);
712  else
713  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
714  }
715 
716  elmat = 0.0;
717  for (int s = 0; s < ir->GetNPoints(); s++)
718  {
719  const IntegrationPoint &ip = ir->IntPoint(s);
720  trial_fe.CalcShape(ip, shape);
721  test_fe.CalcShape(ip, te_shape);
722 
723  Trans.SetIntPoint(&ip);
724  norm = ip.weight * Trans.Weight();
725 
726  MultVWt(te_shape, shape, partelmat);
727 
728  if (VQ)
729  {
730  VQ->Eval(vec, Trans, ip);
731  for (int k = 0; k < vdim; k++)
732  elmat.AddMatrix(norm*vec(k), partelmat, te_nd*k, tr_nd*k);
733  }
734  else if (MQ)
735  {
736  MQ->Eval(mcoeff, Trans, ip);
737  for (int i = 0; i < vdim; i++)
738  for (int j = 0; j < vdim; j++)
739  elmat.AddMatrix(norm*mcoeff(i,j), partelmat, te_nd*i, tr_nd*j);
740  }
741  else
742  {
743  if (Q)
744  norm *= Q->Eval(Trans, ip);
745  partelmat *= norm;
746  for (int k = 0; k < vdim; k++)
747  elmat.AddMatrix(partelmat, te_nd*k, tr_nd*k);
748  }
749  }
750 }
751 
752 void VectorFEDivergenceIntegrator::AssembleElementMatrix2(
753  const FiniteElement &trial_fe, const FiniteElement &test_fe,
754  ElementTransformation &Trans, DenseMatrix &elmat)
755 {
756  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
757 
758 #ifdef MFEM_THREAD_SAFE
759  Vector divshape(trial_nd), shape(test_nd);
760 #else
761  divshape.SetSize(trial_nd);
762  shape.SetSize(test_nd);
763 #endif
764 
765  elmat.SetSize(test_nd, trial_nd);
766 
767  const IntegrationRule *ir = IntRule;
768  if (ir == NULL)
769  {
770  int order = trial_fe.GetOrder() + test_fe.GetOrder() - 1; // <--
771  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
772  }
773 
774  elmat = 0.0;
775  for (i = 0; i < ir->GetNPoints(); i++)
776  {
777  const IntegrationPoint &ip = ir->IntPoint(i);
778  trial_fe.CalcDivShape(ip, divshape);
779  test_fe.CalcShape(ip, shape);
780  double w = ip.weight;
781  if (Q)
782  {
783  Trans.SetIntPoint(&ip);
784  w *= Q->Eval(Trans, ip);
785  }
786  shape *= w;
787  AddMultVWt(shape, divshape, elmat);
788  }
789 }
790 
791 void VectorFECurlIntegrator::AssembleElementMatrix2(
792  const FiniteElement &trial_fe, const FiniteElement &test_fe,
793  ElementTransformation &Trans, DenseMatrix &elmat)
794 {
795  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
796  int dim = trial_fe.GetDim();
797 
798 #ifdef MFEM_THREAD_SAFE
799  DenseMatrix curlshapeTrial(trial_nd, dim);
800  DenseMatrix curlshapeTrial_dFT(trial_nd, dim);
801  DenseMatrix vshapeTest(test_nd, dim);
802 #else
803  curlshapeTrial.SetSize(trial_nd, dim);
804  curlshapeTrial_dFT.SetSize(trial_nd, dim);
805  vshapeTest.SetSize(test_nd, dim);
806 #endif
807 
808  elmat.SetSize(test_nd, trial_nd);
809 
810  const IntegrationRule *ir = IntRule;
811  if (ir == NULL)
812  {
813  int order = trial_fe.GetOrder() + test_fe.GetOrder() - 1; // <--
814  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
815  }
816 
817  elmat = 0.0;
818  for (i = 0; i < ir->GetNPoints(); i++)
819  {
820  const IntegrationPoint &ip = ir->IntPoint(i);
821  Trans.SetIntPoint(&ip);
822  trial_fe.CalcCurlShape(ip, curlshapeTrial);
823  MultABt(curlshapeTrial, Trans.Jacobian(), curlshapeTrial_dFT);
824  test_fe.CalcVShape(Trans, vshapeTest);
825  double w = ip.weight;
826  if (Q)
827  w *= Q->Eval(Trans, ip);
828  vshapeTest *= w;
829  AddMultABt(vshapeTest, curlshapeTrial_dFT, elmat);
830  }
831 }
832 
833 void DerivativeIntegrator::AssembleElementMatrix2 (
834  const FiniteElement &trial_fe,
835  const FiniteElement &test_fe,
836  ElementTransformation &Trans,
837  DenseMatrix &elmat)
838 {
839  int dim = trial_fe.GetDim();
840  int trial_nd = trial_fe.GetDof();
841  int test_nd = test_fe.GetDof();
842 
843  int i, l;
844  double det;
845 
846  elmat.SetSize (test_nd,trial_nd);
847  dshape.SetSize (trial_nd,dim);
848  dshapedxt.SetSize(trial_nd,dim);
849  dshapedxi.SetSize(trial_nd);
850  invdfdx.SetSize(dim);
851  shape.SetSize (test_nd);
852 
853  const IntegrationRule *ir = IntRule;
854  if (ir == NULL)
855  {
856  int order;
857  if (trial_fe.Space() == FunctionSpace::Pk)
858  order = trial_fe.GetOrder() + test_fe.GetOrder() - 1;
859  else
860  order = trial_fe.GetOrder() + test_fe.GetOrder() + dim;
861 
862  if (trial_fe.Space() == FunctionSpace::rQk)
863  ir = &RefinedIntRules.Get(trial_fe.GetGeomType(), order);
864  else
865  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
866  }
867 
868  elmat = 0.0;
869  for(i = 0; i < ir->GetNPoints(); i++) {
870  const IntegrationPoint &ip = ir->IntPoint(i);
871 
872  trial_fe.CalcDShape(ip, dshape);
873 
874  Trans.SetIntPoint (&ip);
875  CalcInverse (Trans.Jacobian(), invdfdx);
876  det = Trans.Weight();
877  Mult (dshape, invdfdx, dshapedxt);
878 
879  test_fe.CalcShape(ip, shape);
880 
881  for (l = 0; l < trial_nd; l++)
882  dshapedxi(l) = dshapedxt(l,xi);
883 
884  shape *= Q.Eval(Trans,ip) * det * ip.weight;
885  AddMultVWt (shape, dshapedxi, elmat);
886  }
887 }
888 
889 void CurlCurlIntegrator::AssembleElementMatrix
891  DenseMatrix &elmat )
892 {
893  int nd = el.GetDof();
894  int dim = el.GetDim();
895  double w;
896 
897 #ifdef MFEM_THREAD_SAFE
898  DenseMatrix Curlshape(nd,dim), Curlshape_dFt(nd,dim);
899 #else
900  Curlshape.SetSize(nd,dim);
901  Curlshape_dFt.SetSize(nd,dim);
902 #endif
903  elmat.SetSize(nd);
904 
905  const IntegrationRule *ir = IntRule;
906  if (ir == NULL)
907  {
908  int order;
909  if (el.Space() == FunctionSpace::Pk)
910  order = 2*el.GetOrder() - 2;
911  else
912  order = 2*el.GetOrder();
913 
914  ir = &IntRules.Get(el.GetGeomType(), order);
915  }
916 
917  elmat = 0.0;
918  for (int i = 0; i < ir->GetNPoints(); i++)
919  {
920  const IntegrationPoint &ip = ir->IntPoint(i);
921  el.CalcCurlShape(ip, Curlshape);
922 
923  Trans.SetIntPoint (&ip);
924 
925  w = ip.weight / Trans.Weight();
926 
927  MultABt(Curlshape, Trans.Jacobian(), Curlshape_dFt);
928 
929  if (Q)
930  w *= Q->Eval(Trans, ip);
931 
932  AddMult_a_AAt(w, Curlshape_dFt, elmat);
933  }
934 }
935 
936 
937 void VectorCurlCurlIntegrator::AssembleElementMatrix(
938  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
939 {
940  int dim = el.GetDim();
941  int dof = el.GetDof();
942  int cld = (dim*(dim-1))/2;
943 
944 #ifdef MFEM_THREAD_SAFE
945  DenseMatrix dshape_hat(dof, dim), dshape(dof, dim);
946  DenseMatrix curlshape(dim*dof, cld), Jadj(dim);
947 #else
948  dshape_hat.SetSize(dof, dim);
949  dshape.SetSize(dof, dim);
950  curlshape.SetSize(dim*dof, cld);
951  Jadj.SetSize(dim);
952 #endif
953 
954  const IntegrationRule *ir = IntRule;
955  if (ir == NULL)
956  {
957  // use the same integration rule as diffusion
958  int order = 2 * Trans.OrderGrad(&el);
959  ir = &IntRules.Get(el.GetGeomType(), order);
960  }
961 
962  elmat = 0.0;
963  for (int i = 0; i < ir->GetNPoints(); i++)
964  {
965  const IntegrationPoint &ip = ir->IntPoint(i);
966  el.CalcDShape(ip, dshape_hat);
967 
968  Trans.SetIntPoint(&ip);
969  CalcAdjugate(Trans.Jacobian(), Jadj);
970  double w = ip.weight / Trans.Weight();
971 
972  Mult(dshape_hat, Jadj, dshape);
973  dshape.GradToCurl(curlshape);
974 
975  if (Q)
976  w *= Q->Eval(Trans, ip);
977 
978  AddMult_a_AAt(w, curlshape, elmat);
979  }
980 }
981 
982 double VectorCurlCurlIntegrator::GetElementEnergy(
983  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
984 {
985  int dim = el.GetDim();
986  int dof = el.GetDof();
987 
988 #ifdef MFEM_THREAD_SAFE
989  DenseMatrix dshape_hat(dof, dim), Jadj(dim), grad_hat(dim), grad(dim);
990 #else
991  dshape_hat.SetSize(dof, dim);
992  Jadj.SetSize(dim);
993  grad_hat.SetSize(dim);
994  grad.SetSize(dim);
995 #endif
996  DenseMatrix elfun_mat(elfun.GetData(), dof, dim);
997 
998  const IntegrationRule *ir = IntRule;
999  if (ir == NULL)
1000  {
1001  // use the same integration rule as diffusion
1002  int order = 2 * Tr.OrderGrad(&el);
1003  ir = &IntRules.Get(el.GetGeomType(), order);
1004  }
1005 
1006  double energy = 0.;
1007  for (int i = 0; i < ir->GetNPoints(); i++)
1008  {
1009  const IntegrationPoint &ip = ir->IntPoint(i);
1010  el.CalcDShape(ip, dshape_hat);
1011 
1012  MultAtB(elfun_mat, dshape_hat, grad_hat);
1013 
1014  Tr.SetIntPoint(&ip);
1015  CalcAdjugate(Tr.Jacobian(), Jadj);
1016  double w = ip.weight / Tr.Weight();
1017 
1018  Mult(grad_hat, Jadj, grad);
1019 
1020  if (dim == 2)
1021  {
1022  double curl = grad(0,1) - grad(1,0);
1023  w *= curl * curl;
1024  }
1025  else
1026  {
1027  double curl_x = grad(2,1) - grad(1,2);
1028  double curl_y = grad(0,2) - grad(2,0);
1029  double curl_z = grad(1,0) - grad(0,1);
1030  w *= curl_x * curl_x + curl_y * curl_y + curl_z * curl_z;
1031  }
1032 
1033  if (Q)
1034  w *= Q->Eval(Tr, ip);
1035 
1036  energy += w;
1037  }
1038 
1039  elfun_mat.ClearExternalData();
1040 
1041  return 0.5 * energy;
1042 }
1043 
1044 
1045 void VectorFEMassIntegrator::AssembleElementMatrix(
1046  const FiniteElement &el,
1047  ElementTransformation &Trans,
1048  DenseMatrix &elmat)
1049 {
1050  int dof = el.GetDof();
1051  int dim = el.GetDim();
1052 
1053  double w;
1054 
1055 #ifdef MFEM_THREAD_SAFE
1056  Vector D(VQ ? VQ->GetVDim() : 0);
1057  DenseMatrix vshape(dof, dim);
1058  DenseMatrix K(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
1059 #else
1060  vshape.SetSize(dof,dim);
1061  D.SetSize(VQ ? VQ->GetVDim() : 0);
1062  K.SetSize(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
1063 #endif
1064  DenseMatrix tmp(vshape.Height(), K.Width());
1065 
1066  elmat.SetSize(dof);
1067  elmat = 0.0;
1068 
1069  const IntegrationRule *ir = IntRule;
1070  if (ir == NULL)
1071  {
1072  // int order = 2 * el.GetOrder();
1073  int order = Trans.OrderW() + 2 * el.GetOrder();
1074  ir = &IntRules.Get(el.GetGeomType(), order);
1075  }
1076 
1077  for (int i = 0; i < ir->GetNPoints(); i++)
1078  {
1079  const IntegrationPoint &ip = ir->IntPoint(i);
1080 
1081  Trans.SetIntPoint (&ip);
1082 
1083  el.CalcVShape(Trans, vshape);
1084 
1085  w = ip.weight * Trans.Weight();
1086  if (MQ)
1087  {
1088  MQ->Eval(K, Trans, ip);
1089  K *= w;
1090  Mult(vshape,K,tmp);
1091  AddMultABt(tmp,vshape,elmat);
1092  }
1093  else if (VQ)
1094  {
1095  VQ->Eval(D, Trans, ip);
1096  D *= w;
1097  AddMultADAt(vshape, D, elmat);
1098  }
1099  else
1100  {
1101  if (Q)
1102  w *= Q -> Eval (Trans, ip);
1103  AddMult_a_AAt (w, vshape, elmat);
1104  }
1105  }
1106 }
1107 
1108 void VectorFEMassIntegrator::AssembleElementMatrix2(
1109  const FiniteElement &trial_fe, const FiniteElement &test_fe,
1110  ElementTransformation &Trans, DenseMatrix &elmat)
1111 {
1112  // assume test_fe is scalar FE and trial_fe is vector FE
1113  int dim = test_fe.GetDim();
1114  int trial_dof = trial_fe.GetDof();
1115  int test_dof = test_fe.GetDof();
1116  double w;
1117 
1118  if (VQ || MQ)
1119  mfem_error("VectorFEMassIntegrator::AssembleElementMatrix2(...)\n"
1120  " is not implemented for vector/tensor permeability");
1121 
1122 #ifdef MFEM_THREAD_SAFE
1123  DenseMatrix vshape(trial_dof, dim);
1124  Vector shape(test_dof);
1125 #else
1126  vshape.SetSize(trial_dof, dim);
1127  shape.SetSize(test_dof);
1128 #endif
1129 
1130  elmat.SetSize (dim*test_dof, trial_dof);
1131 
1132  const IntegrationRule *ir = IntRule;
1133  if (ir == NULL)
1134  {
1135  int order = (Trans.OrderW() + test_fe.GetOrder() + trial_fe.GetOrder());
1136  ir = &IntRules.Get(test_fe.GetGeomType(), order);
1137  }
1138 
1139  elmat = 0.0;
1140  for (int i = 0; i < ir->GetNPoints(); i++)
1141  {
1142  const IntegrationPoint &ip = ir->IntPoint(i);
1143 
1144  Trans.SetIntPoint (&ip);
1145 
1146  trial_fe.CalcVShape(Trans, vshape);
1147  test_fe.CalcShape(ip, shape);
1148 
1149  w = ip.weight * Trans.Weight();
1150  if (Q)
1151  w *= Q -> Eval (Trans, ip);
1152 
1153  for (int d = 0; d < dim; d++)
1154  {
1155  for (int j = 0; j < test_dof; j++)
1156  {
1157  for (int k = 0; k < trial_dof; k++)
1158  {
1159  elmat(d * test_dof + j, k) += w * shape(j) * vshape(k, d);
1160  }
1161  }
1162  }
1163  }
1164 }
1165 
1166 
1167 void VectorDivergenceIntegrator::AssembleElementMatrix2(
1168  const FiniteElement &trial_fe,
1169  const FiniteElement &test_fe,
1170  ElementTransformation &Trans,
1171  DenseMatrix &elmat)
1172 {
1173  int dim = trial_fe.GetDim();
1174  int trial_dof = trial_fe.GetDof();
1175  int test_dof = test_fe.GetDof();
1176  double c;
1177 
1178  dshape.SetSize (trial_dof, dim);
1179  gshape.SetSize (trial_dof, dim);
1180  Jadj.SetSize (dim);
1181  divshape.SetSize (dim*trial_dof);
1182  shape.SetSize (test_dof);
1183 
1184  elmat.SetSize (test_dof, dim*trial_dof);
1185 
1186  const IntegrationRule *ir = IntRule;
1187  if (ir == NULL)
1188  {
1189  int order = Trans.OrderGrad(&trial_fe) + test_fe.GetOrder();
1190  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
1191  }
1192 
1193  elmat = 0.0;
1194 
1195  for (int i = 0; i < ir -> GetNPoints(); i++)
1196  {
1197  const IntegrationPoint &ip = ir->IntPoint(i);
1198 
1199  trial_fe.CalcDShape (ip, dshape);
1200  test_fe.CalcShape (ip, shape);
1201 
1202  Trans.SetIntPoint (&ip);
1203  CalcAdjugate(Trans.Jacobian(), Jadj);
1204 
1205  Mult (dshape, Jadj, gshape);
1206 
1207  gshape.GradToDiv (divshape);
1208 
1209  c = ip.weight;
1210  if (Q)
1211  c *= Q -> Eval (Trans, ip);
1212 
1213  // elmat += c * shape * divshape ^ t
1214  shape *= c;
1215  AddMultVWt (shape, divshape, elmat);
1216  }
1217 }
1218 
1219 
1220 void DivDivIntegrator::AssembleElementMatrix(
1221  const FiniteElement &el,
1222  ElementTransformation &Trans,
1223  DenseMatrix &elmat)
1224 {
1225  int dof = el.GetDof();
1226  double c;
1227 
1228 #ifdef MFEM_THREAD_SAFE
1229  Vector divshape(dof);
1230 #else
1231  divshape.SetSize(dof);
1232 #endif
1233  elmat.SetSize(dof);
1234 
1235  const IntegrationRule *ir = IntRule;
1236  if (ir == NULL)
1237  {
1238  int order = 2 * el.GetOrder() - 2; // <--- OK for RTk
1239  ir = &IntRules.Get(el.GetGeomType(), order);
1240  }
1241 
1242  elmat = 0.0;
1243 
1244  for (int i = 0; i < ir -> GetNPoints(); i++)
1245  {
1246  const IntegrationPoint &ip = ir->IntPoint(i);
1247 
1248  el.CalcDivShape (ip, divshape);
1249 
1250  Trans.SetIntPoint (&ip);
1251  c = ip.weight / Trans.Weight();
1252 
1253  if (Q)
1254  c *= Q -> Eval (Trans, ip);
1255 
1256  // elmat += c * divshape * divshape ^ t
1257  AddMult_a_VVt (c, divshape, elmat);
1258  }
1259 }
1260 
1261 
1262 void VectorDiffusionIntegrator::AssembleElementMatrix(
1263  const FiniteElement &el,
1264  ElementTransformation &Trans,
1265  DenseMatrix &elmat)
1266 {
1267  int dim = el.GetDim();
1268  int dof = el.GetDof();
1269 
1270  double norm;
1271 
1272  elmat.SetSize (dim * dof);
1273 
1274  Jinv. SetSize (dim);
1275  dshape.SetSize (dof, dim);
1276  gshape.SetSize (dof, dim);
1277  pelmat.SetSize (dof);
1278 
1279  const IntegrationRule *ir = IntRule;
1280  if (ir == NULL)
1281  {
1282  // integrant is rational function if det(J) is not constant
1283  int order = 2 * Trans.OrderGrad(&el); // order of the numerator
1284  if (el.Space() == FunctionSpace::rQk)
1285  ir = &RefinedIntRules.Get(el.GetGeomType(), order);
1286  else
1287  ir = &IntRules.Get(el.GetGeomType(), order);
1288  }
1289 
1290  elmat = 0.0;
1291 
1292  for (int i = 0; i < ir -> GetNPoints(); i++)
1293  {
1294  const IntegrationPoint &ip = ir->IntPoint(i);
1295 
1296  el.CalcDShape (ip, dshape);
1297 
1298  Trans.SetIntPoint (&ip);
1299  norm = ip.weight * Trans.Weight();
1300  CalcInverse (Trans.Jacobian(), Jinv);
1301 
1302  Mult (dshape, Jinv, gshape);
1303 
1304  MultAAt (gshape, pelmat);
1305 
1306  if (Q)
1307  norm *= Q -> Eval (Trans, ip);
1308 
1309  pelmat *= norm;
1310 
1311  for (int d = 0; d < dim; d++)
1312  {
1313  for (int k = 0; k < dof; k++)
1314  for (int l = 0; l < dof; l++)
1315  elmat (dof*d+k, dof*d+l) += pelmat (k, l);
1316  }
1317  }
1318 }
1319 
1320 
1321 void ElasticityIntegrator::AssembleElementMatrix(
1322  const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
1323 {
1324  int dof = el.GetDof();
1325  int dim = el.GetDim();
1326  double w, L, M;
1327 
1328 #ifdef MFEM_THREAD_SAFE
1329  DenseMatrix dshape(dof, dim), Jinv(dim), gshape(dof, dim), pelmat(dof);
1330  Vector divshape(dim*dof);
1331 #else
1332  Jinv.SetSize(dim);
1333  dshape.SetSize(dof, dim);
1334  gshape.SetSize(dof, dim);
1335  pelmat.SetSize(dof);
1336  divshape.SetSize(dim*dof);
1337 #endif
1338 
1339  elmat.SetSize(dof * dim);
1340 
1341  const IntegrationRule *ir = IntRule;
1342  if (ir == NULL)
1343  {
1344  int order = 2 * Trans.OrderGrad(&el); // correct order?
1345  ir = &IntRules.Get(el.GetGeomType(), order);
1346  }
1347 
1348  elmat = 0.0;
1349 
1350  for (int i = 0; i < ir -> GetNPoints(); i++)
1351  {
1352  const IntegrationPoint &ip = ir->IntPoint(i);
1353 
1354  el.CalcDShape(ip, dshape);
1355 
1356  Trans.SetIntPoint(&ip);
1357  w = ip.weight * Trans.Weight();
1358  CalcInverse(Trans.Jacobian(), Jinv);
1359  Mult(dshape, Jinv, gshape);
1360  MultAAt(gshape, pelmat);
1361  gshape.GradToDiv (divshape);
1362 
1363  M = mu->Eval(Trans, ip);
1364  if (lambda)
1365  L = lambda->Eval(Trans, ip);
1366  else
1367  {
1368  L = q_lambda * M;
1369  M = q_mu * M;
1370  }
1371 
1372  if (L != 0.0)
1373  AddMult_a_VVt(L * w, divshape, elmat);
1374 
1375  if (M != 0.0)
1376  {
1377  for (int d = 0; d < dim; d++)
1378  {
1379  for (int k = 0; k < dof; k++)
1380  for (int l = 0; l < dof; l++)
1381  elmat (dof*d+k, dof*d+l) += (M * w) * pelmat(k, l);
1382  }
1383  for (int i = 0; i < dim; i++)
1384  for (int j = 0; j < dim; j++)
1385  {
1386  for (int k = 0; k < dof; k++)
1387  for (int l = 0; l < dof; l++)
1388  elmat(dof*i+k, dof*j+l) +=
1389  (M * w) * gshape(k, j) * gshape(l, i);
1390  // + (L * w) * gshape(k, i) * gshape(l, j)
1391  }
1392  }
1393  }
1394 }
1395 
1396 void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
1397  const FiniteElement &el2,
1399  DenseMatrix &elmat)
1400 {
1401  int dim, ndof1, ndof2;
1402 
1403  double un, a, b, w;
1404 
1405  dim = el1.GetDim();
1406  ndof1 = el1.GetDof();
1407  Vector vu(dim), nor(dim);
1408 
1409  if (Trans.Elem2No >= 0)
1410  ndof2 = el2.GetDof();
1411  else
1412  ndof2 = 0;
1413 
1414  shape1.SetSize(ndof1);
1415  shape2.SetSize(ndof2);
1416  elmat.SetSize(ndof1 + ndof2);
1417  elmat = 0.0;
1418 
1419  const IntegrationRule *ir = IntRule;
1420  if (ir == NULL)
1421  {
1422  int order;
1423  // Assuming order(u)==order(mesh)
1424  if (Trans.Elem2No >= 0)
1425  order = (min(Trans.Elem1->OrderW(), Trans.Elem2->OrderW()) +
1426  2*max(el1.GetOrder(), el2.GetOrder()));
1427  else
1428  order = Trans.Elem1->OrderW() + 2*el1.GetOrder();
1429  if (el1.Space() == FunctionSpace::Pk)
1430  order++;
1431  ir = &IntRules.Get(Trans.FaceGeom, order);
1432  }
1433 
1434  for (int p = 0; p < ir->GetNPoints(); p++)
1435  {
1436  const IntegrationPoint &ip = ir->IntPoint(p);
1437  IntegrationPoint eip1, eip2;
1438  Trans.Loc1.Transform(ip, eip1);
1439  if (ndof2)
1440  Trans.Loc2.Transform(ip, eip2);
1441  el1.CalcShape(eip1, shape1);
1442 
1443  Trans.Face->SetIntPoint(&ip);
1444  Trans.Elem1->SetIntPoint(&eip1);
1445 
1446  u->Eval(vu, *Trans.Elem1, eip1);
1447 
1448  if (dim == 1)
1449  nor(0) = 2*eip1.x - 1.0;
1450  else
1451  CalcOrtho(Trans.Face->Jacobian(), nor);
1452 
1453  un = vu * nor;
1454  a = 0.5 * alpha * un;
1455  b = beta * fabs(un);
1456  // note: if |alpha/2|==|beta| then |a|==|b|, i.e. (a==b) or (a==-b)
1457  // and therefore two blocks in the element matrix contribution
1458  // (from the current quadrature point) are 0
1459 
1460  if (rho)
1461  {
1462  double rho_p;
1463  if (un >= 0.0 && ndof2)
1464  {
1465  Trans.Elem2->SetIntPoint(&eip2);
1466  rho_p = rho->Eval(*Trans.Elem2, eip2);
1467  }
1468  else
1469  {
1470  rho_p = rho->Eval(*Trans.Elem1, eip1);
1471  }
1472  a *= rho_p;
1473  b *= rho_p;
1474  }
1475 
1476  w = ip.weight * (a+b);
1477  if (w != 0.0)
1478  {
1479  for (int i = 0; i < ndof1; i++)
1480  for (int j = 0; j < ndof1; j++)
1481  elmat(i, j) += w * shape1(i) * shape1(j);
1482  }
1483 
1484  if (ndof2)
1485  {
1486  el2.CalcShape(eip2, shape2);
1487 
1488  if (w != 0.0)
1489  for (int i = 0; i < ndof2; i++)
1490  for (int j = 0; j < ndof1; j++)
1491  elmat(ndof1+i, j) -= w * shape2(i) * shape1(j);
1492 
1493  w = ip.weight * (b-a);
1494  if (w != 0.0)
1495  {
1496  for (int i = 0; i < ndof2; i++)
1497  for (int j = 0; j < ndof2; j++)
1498  elmat(ndof1+i, ndof1+j) += w * shape2(i) * shape2(j);
1499 
1500  for (int i = 0; i < ndof1; i++)
1501  for (int j = 0; j < ndof2; j++)
1502  elmat(i, ndof1+j) -= w * shape1(i) * shape2(j);
1503  }
1504  }
1505  }
1506 }
1507 
1508 void DGDiffusionIntegrator::AssembleFaceMatrix(
1509  const FiniteElement &el1, const FiniteElement &el2,
1510  FaceElementTransformations &Trans, DenseMatrix &elmat)
1511 {
1512  int dim, ndof1, ndof2, ndofs;
1513  bool kappa_is_nonzero = (kappa != 0.);
1514  double w, wq = 0.0;
1515 
1516  dim = el1.GetDim();
1517  ndof1 = el1.GetDof();
1518 
1519  nor.SetSize(dim);
1520  nh.SetSize(dim);
1521  ni.SetSize(dim);
1522  adjJ.SetSize(dim);
1523  if (MQ)
1524  mq.SetSize(dim);
1525 
1526  shape1.SetSize(ndof1);
1527  dshape1.SetSize(ndof1, dim);
1528  dshape1dn.SetSize(ndof1);
1529  if (Trans.Elem2No >= 0)
1530  {
1531  ndof2 = el2.GetDof();
1532  shape2.SetSize(ndof2);
1533  dshape2.SetSize(ndof2, dim);
1534  dshape2dn.SetSize(ndof2);
1535  }
1536  else
1537  ndof2 = 0;
1538 
1539  ndofs = ndof1 + ndof2;
1540  elmat.SetSize(ndofs);
1541  elmat = 0.0;
1542  if (kappa_is_nonzero)
1543  {
1544  jmat.SetSize(ndofs);
1545  jmat = 0.;
1546  }
1547 
1548  const IntegrationRule *ir = IntRule;
1549  if (ir == NULL)
1550  {
1551  // a simple choice for the integration order; is this OK?
1552  int order;
1553  if (ndof2)
1554  order = 2*max(el1.GetOrder(), el2.GetOrder());
1555  else
1556  order = 2*el1.GetOrder();
1557  ir = &IntRules.Get(Trans.FaceGeom, order);
1558  }
1559 
1560  // assmble: < {(Q \nabla u).n},[v] > --> elmat
1561  // kappa < {h^{-1} Q} [u],[v] > --> jmat
1562  for (int p = 0; p < ir->GetNPoints(); p++)
1563  {
1564  const IntegrationPoint &ip = ir->IntPoint(p);
1565  IntegrationPoint eip1, eip2;
1566 
1567  Trans.Loc1.Transform(ip, eip1);
1568  Trans.Face->SetIntPoint(&ip);
1569  if (dim == 1)
1570  nor(0) = 2*eip1.x - 1.0;
1571  else
1572  CalcOrtho(Trans.Face->Jacobian(), nor);
1573 
1574  el1.CalcShape(eip1, shape1);
1575  el1.CalcDShape(eip1, dshape1);
1576  Trans.Elem1->SetIntPoint(&eip1);
1577  w = ip.weight/Trans.Elem1->Weight();
1578  if (ndof2)
1579  w /= 2;
1580  if (!MQ)
1581  {
1582  if (Q)
1583  w *= Q->Eval(*Trans.Elem1, eip1);
1584  ni.Set(w, nor);
1585  }
1586  else
1587  {
1588  nh.Set(w, nor);
1589  MQ->Eval(mq, *Trans.Elem1, eip1);
1590  mq.MultTranspose(nh, ni);
1591  }
1592  CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
1593  adjJ.Mult(ni, nh);
1594  if (kappa_is_nonzero)
1595  wq = ni * nor;
1596  // Note: in the jump term, we use 1/h1 = |nor|/det(J1) which is
1597  // independent of Loc1 and always gives the size of element 1 in
1598  // direction perpendicular to the face. Indeed, for linear transformation
1599  // |nor|=measure(face)/measure(ref. face),
1600  // det(J1)=measure(element)/measure(ref. element),
1601  // and the ratios measure(ref. element)/measure(ref. face) are
1602  // compatible for all element/face pairs.
1603  // For example: meas(ref. tetrahedron)/meas(ref. triangle) = 1/3, and
1604  // for any tetrahedron vol(tet)=(1/3)*height*area(base).
1605  // For interior faces: q_e/h_e=(q1/h1+q2/h2)/2.
1606 
1607  dshape1.Mult(nh, dshape1dn);
1608  for (int i = 0; i < ndof1; i++)
1609  for (int j = 0; j < ndof1; j++)
1610  elmat(i, j) += shape1(i) * dshape1dn(j);
1611 
1612  if (ndof2)
1613  {
1614  Trans.Loc2.Transform(ip, eip2);
1615  el2.CalcShape(eip2, shape2);
1616  el2.CalcDShape(eip2, dshape2);
1617  Trans.Elem2->SetIntPoint(&eip2);
1618  w = ip.weight/2/Trans.Elem2->Weight();
1619  if (!MQ)
1620  {
1621  if (Q)
1622  w *= Q->Eval(*Trans.Elem2, eip2);
1623  ni.Set(w, nor);
1624  }
1625  else
1626  {
1627  nh.Set(w, nor);
1628  MQ->Eval(mq, *Trans.Elem2, eip2);
1629  mq.MultTranspose(nh, ni);
1630  }
1631  CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
1632  adjJ.Mult(ni, nh);
1633  if (kappa_is_nonzero)
1634  wq += ni * nor;
1635 
1636  dshape2.Mult(nh, dshape2dn);
1637 
1638  for (int i = 0; i < ndof1; i++)
1639  for (int j = 0; j < ndof2; j++)
1640  elmat(i, ndof1 + j) += shape1(i) * dshape2dn(j);
1641 
1642  for (int i = 0; i < ndof2; i++)
1643  for (int j = 0; j < ndof1; j++)
1644  elmat(ndof1 + i, j) -= shape2(i) * dshape1dn(j);
1645 
1646  for (int i = 0; i < ndof2; i++)
1647  for (int j = 0; j < ndof2; j++)
1648  elmat(ndof1 + i, ndof1 + j) -= shape2(i) * dshape2dn(j);
1649  }
1650 
1651  if (kappa_is_nonzero)
1652  {
1653  // only assemble the lower triangular part of jmat
1654  wq *= kappa;
1655  for (int i = 0; i < ndof1; i++)
1656  {
1657  const double wsi = wq*shape1(i);
1658  for (int j = 0; j <= i; j++)
1659  jmat(i, j) += wsi * shape1(j);
1660  }
1661  if (ndof2)
1662  {
1663  for (int i = 0; i < ndof2; i++)
1664  {
1665  const int i2 = ndof1 + i;
1666  const double wsi = wq*shape2(i);
1667  for (int j = 0; j < ndof1; j++)
1668  jmat(i2, j) -= wsi * shape1(j);
1669  for (int j = 0; j <= i; j++)
1670  jmat(i2, ndof1 + j) += wsi * shape2(j);
1671  }
1672  }
1673  }
1674  }
1675 
1676  // elmat := -elmat + sigma*elmat^t + jmat
1677  if (kappa_is_nonzero)
1678  {
1679  for (int i = 0; i < ndofs; i++)
1680  {
1681  for (int j = 0; j < i; j++)
1682  {
1683  double aij = elmat(i,j), aji = elmat(j,i), mij = jmat(i,j);
1684  elmat(i,j) = sigma*aji - aij + mij;
1685  elmat(j,i) = sigma*aij - aji + mij;
1686  }
1687  elmat(i,i) = (sigma - 1.)*elmat(i,i) + jmat(i,i);
1688  }
1689  }
1690  else
1691  {
1692  for (int i = 0; i < ndofs; i++)
1693  {
1694  for (int j = 0; j < i; j++)
1695  {
1696  double aij = elmat(i,j), aji = elmat(j,i);
1697  elmat(i,j) = sigma*aji - aij;
1698  elmat(j,i) = sigma*aij - aji;
1699  }
1700  elmat(i,i) *= (sigma - 1.);
1701  }
1702  }
1703 }
1704 
1705 void TraceJumpIntegrator::AssembleFaceMatrix(
1706  const FiniteElement &trial_face_fe, const FiniteElement &test_fe1,
1707  const FiniteElement &test_fe2, FaceElementTransformations &Trans,
1708  DenseMatrix &elmat)
1709 {
1710  int i, j, face_ndof, ndof1, ndof2;
1711  int order;
1712 
1713  double w;
1714 
1715  face_ndof = trial_face_fe.GetDof();
1716  ndof1 = test_fe1.GetDof();
1717 
1718  face_shape.SetSize(face_ndof);
1719  shape1.SetSize(ndof1);
1720 
1721  if (Trans.Elem2No >= 0)
1722  {
1723  ndof2 = test_fe2.GetDof();
1724  shape2.SetSize(ndof2);
1725  }
1726  else
1727  ndof2 = 0;
1728 
1729  elmat.SetSize(ndof1 + ndof2, face_ndof);
1730  elmat = 0.0;
1731 
1732  const IntegrationRule *ir = IntRule;
1733  if (ir == NULL)
1734  {
1735  if (Trans.Elem2No >= 0)
1736  order = max(test_fe1.GetOrder(), test_fe2.GetOrder());
1737  else
1738  order = test_fe1.GetOrder();
1739  order += trial_face_fe.GetOrder();
1740  if (trial_face_fe.GetMapType() == FiniteElement::VALUE)
1741  order += Trans.Face->OrderW();
1742  ir = &IntRules.Get(Trans.FaceGeom, order);
1743  }
1744 
1745  for (int p = 0; p < ir->GetNPoints(); p++)
1746  {
1747  const IntegrationPoint &ip = ir->IntPoint(p);
1748  IntegrationPoint eip1, eip2;
1749  // Trace finite element shape function
1750  Trans.Face->SetIntPoint(&ip);
1751  trial_face_fe.CalcShape(ip, face_shape);
1752  // Side 1 finite element shape function
1753  Trans.Loc1.Transform(ip, eip1);
1754  test_fe1.CalcShape(eip1, shape1);
1755  Trans.Elem1->SetIntPoint(&eip1);
1756  if (ndof2)
1757  {
1758  // Side 2 finite element shape function
1759  Trans.Loc2.Transform(ip, eip2);
1760  test_fe2.CalcShape(eip2, shape2);
1761  Trans.Elem2->SetIntPoint(&eip2);
1762  }
1763  w = ip.weight;
1764  if (trial_face_fe.GetMapType() == FiniteElement::VALUE)
1765  w *= Trans.Face->Weight();
1766  face_shape *= w;
1767  for (i = 0; i < ndof1; i++)
1768  for (j = 0; j < face_ndof; j++)
1769  elmat(i, j) += shape1(i) * face_shape(j);
1770  if (ndof2)
1771  {
1772  // Subtract contribution from side 2
1773  for (i = 0; i < ndof2; i++)
1774  for (j = 0; j < face_ndof; j++)
1775  elmat(ndof1+i, j) -= shape2(i) * face_shape(j);
1776  }
1777  }
1778 }
1779 
1780 }
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
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
Definition: densemat.cpp:2782
int GetDim() const
Returns the space dimension for the finite element.
Definition: fe.hpp:80
ElementTransformation * Face
Definition: eltrans.hpp:101
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3025
Class for integration rule.
Definition: intrules.hpp:63
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Definition: fe.cpp:30
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:3007
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:194
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:248
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:351
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:202
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
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:102
int Space() const
Returns the type of space on each element.
Definition: fe.hpp:92
int GetMapType() const
Definition: fe.hpp:96
ElementTransformation * Elem2
Definition: eltrans.hpp:101
double * GetData() const
Definition: vector.hpp:80
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2710
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 CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:465
const IntegrationRule & GetNodes() const
Definition: fe.hpp:110
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:3059
virtual double Weight()=0
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2588
int GetGeomType() const
Returns the geometry type:
Definition: fe.hpp:83
const double kappa
Definition: ex3.cpp:186
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:2036
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:2998
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2859
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
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2732
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Definition: fe.cpp:44
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0&lt;=i&lt;A.Height, 0&lt;=j&lt;A.Width.
Definition: densemat.cpp:2246
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Definition: fe.cpp:51
Class for integration point with weight.
Definition: intrules.hpp:25
virtual int OrderGrad(const FiniteElement *fe)=0
order of adj(J)^t.grad(fi)
void GradToCurl(DenseMatrix &curl)
Definition: densemat.cpp:2096
ElementTransformation * Elem1
Definition: eltrans.hpp:101
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2912
Vector data type.
Definition: vector.hpp:29
IntegrationRules RefinedIntRules(1)
A global object with all refined integration rules.
Definition: intrules.hpp:267
void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt += A D A^t, where D is diagonal.
Definition: densemat.cpp:2744
virtual int GetSpaceDim()=0
virtual const DenseMatrix & Jacobian()=0
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
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:2965