MFEM  v4.5.2
Finite element discretization library
bilininteg.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 // 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::AssemblePA(const FiniteElementSpace&)
24 {
25  mfem_error ("BilinearFormIntegrator::AssemblePA(fes)\n"
26  " is not implemented for this class.");
27 }
28 
29 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&,
30  const FiniteElementSpace&)
31 {
32  mfem_error ("BilinearFormIntegrator::AssemblePA(fes, fes)\n"
33  " is not implemented for this class.");
34 }
35 
36 void BilinearFormIntegrator::AssemblePAInteriorFaces(const FiniteElementSpace&)
37 {
38  mfem_error ("BilinearFormIntegrator::AssemblePAInteriorFaces(...)\n"
39  " is not implemented for this class.");
40 }
41 
42 void BilinearFormIntegrator::AssemblePABoundaryFaces(const FiniteElementSpace&)
43 {
44  mfem_error ("BilinearFormIntegrator::AssemblePABoundaryFaces(...)\n"
45  " is not implemented for this class.");
46 }
47 
48 void BilinearFormIntegrator::AssembleDiagonalPA(Vector &)
49 {
50  mfem_error ("BilinearFormIntegrator::AssembleDiagonalPA(...)\n"
51  " is not implemented for this class.");
52 }
53 
54 void BilinearFormIntegrator::AssembleEA(const FiniteElementSpace &fes,
55  Vector &emat,
56  const bool add)
57 {
58  mfem_error ("BilinearFormIntegrator::AssembleEA(...)\n"
59  " is not implemented for this class.");
60 }
61 
62 void BilinearFormIntegrator::AssembleEAInteriorFaces(const FiniteElementSpace
63  &fes,
64  Vector &ea_data_int,
65  Vector &ea_data_ext,
66  const bool add)
67 {
68  mfem_error ("BilinearFormIntegrator::AssembleEAInteriorFaces(...)\n"
69  " is not implemented for this class.");
70 }
71 
72 void BilinearFormIntegrator::AssembleEABoundaryFaces(const FiniteElementSpace
73  &fes,
74  Vector &ea_data_bdr,
75  const bool add)
76 {
77  mfem_error ("BilinearFormIntegrator::AssembleEABoundaryFaces(...)\n"
78  " is not implemented for this class.");
79 }
80 
81 void BilinearFormIntegrator::AssembleDiagonalPA_ADAt(const Vector &, Vector &)
82 {
83  MFEM_ABORT("BilinearFormIntegrator::AssembleDiagonalPA_ADAt(...)\n"
84  " is not implemented for this class.");
85 }
86 
87 void BilinearFormIntegrator::AddMultPA(const Vector &, Vector &) const
88 {
89  mfem_error ("BilinearFormIntegrator::MultAssembled(...)\n"
90  " is not implemented for this class.");
91 }
92 
93 void BilinearFormIntegrator::AddMultTransposePA(const Vector &, Vector &) const
94 {
95  mfem_error ("BilinearFormIntegrator::AddMultTransposePA(...)\n"
96  " is not implemented for this class.");
97 }
98 
99 void BilinearFormIntegrator::AssembleMF(const FiniteElementSpace &fes)
100 {
101  mfem_error ("BilinearFormIntegrator::AssembleMF(...)\n"
102  " is not implemented for this class.");
103 }
104 
105 void BilinearFormIntegrator::AddMultMF(const Vector &, Vector &) const
106 {
107  mfem_error ("BilinearFormIntegrator::AddMultMF(...)\n"
108  " is not implemented for this class.");
109 }
110 
111 void BilinearFormIntegrator::AddMultTransposeMF(const Vector &, Vector &) const
112 {
113  mfem_error ("BilinearFormIntegrator::AddMultTransposeMF(...)\n"
114  " is not implemented for this class.");
115 }
116 
117 void BilinearFormIntegrator::AssembleDiagonalMF(Vector &)
118 {
119  mfem_error ("BilinearFormIntegrator::AssembleDiagonalMF(...)\n"
120  " is not implemented for this class.");
121 }
122 
123 void BilinearFormIntegrator::AssembleElementMatrix (
125  DenseMatrix &elmat )
126 {
127  mfem_error ("BilinearFormIntegrator::AssembleElementMatrix(...)\n"
128  " is not implemented for this class.");
129 }
130 
131 void BilinearFormIntegrator::AssembleElementMatrix2 (
132  const FiniteElement &el1, const FiniteElement &el2,
134 {
135  mfem_error ("BilinearFormIntegrator::AssembleElementMatrix2(...)\n"
136  " is not implemented for this class.");
137 }
138 
139 void BilinearFormIntegrator::AssembleFaceMatrix (
140  const FiniteElement &el1, const FiniteElement &el2,
142 {
143  mfem_error ("BilinearFormIntegrator::AssembleFaceMatrix(...)\n"
144  " is not implemented for this class.");
145 }
146 
147 void BilinearFormIntegrator::AssembleFaceMatrix(
148  const FiniteElement &trial_face_fe, const FiniteElement &test_fe1,
150  DenseMatrix &elmat)
151 {
152  MFEM_ABORT("AssembleFaceMatrix (mixed form) is not implemented for this"
153  " Integrator class.");
154 }
155 
156 void BilinearFormIntegrator::AssembleElementVector(
157  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
158  Vector &elvect)
159 {
160  // Note: This default implementation is general but not efficient
161  DenseMatrix elmat;
162  AssembleElementMatrix(el, Tr, elmat);
163  elvect.SetSize(elmat.Height());
164  elmat.Mult(elfun, elvect);
165 }
166 
167 void BilinearFormIntegrator::AssembleFaceVector(
168  const FiniteElement &el1, const FiniteElement &el2,
169  FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect)
170 {
171  // Note: This default implementation is general but not efficient
172  DenseMatrix elmat;
173  AssembleFaceMatrix(el1, el2, Tr, elmat);
174  elvect.SetSize(elmat.Height());
175  elmat.Mult(elfun, elvect);
176 }
177 
178 void TransposeIntegrator::SetIntRule(const IntegrationRule *ir)
179 {
180  IntRule = ir;
181  bfi->SetIntRule(ir);
182 }
183 
184 void TransposeIntegrator::AssembleElementMatrix (
186 {
187  bfi -> AssembleElementMatrix (el, Trans, bfi_elmat);
188  // elmat = bfi_elmat^t
189  elmat.Transpose (bfi_elmat);
190 }
191 
192 void TransposeIntegrator::AssembleElementMatrix2 (
193  const FiniteElement &trial_fe, const FiniteElement &test_fe,
195 {
196  bfi -> AssembleElementMatrix2 (test_fe, trial_fe, Trans, bfi_elmat);
197  // elmat = bfi_elmat^t
198  elmat.Transpose (bfi_elmat);
199 }
200 
201 void TransposeIntegrator::AssembleFaceMatrix (
202  const FiniteElement &el1, const FiniteElement &el2,
204 {
205  bfi -> AssembleFaceMatrix (el1, el2, Trans, bfi_elmat);
206  // elmat = bfi_elmat^t
207  elmat.Transpose (bfi_elmat);
208 }
209 
210 void LumpedIntegrator::SetIntRule(const IntegrationRule *ir)
211 {
212  IntRule = ir;
213  bfi->SetIntRule(ir);
214 }
215 
216 void LumpedIntegrator::AssembleElementMatrix (
218 {
219  bfi -> AssembleElementMatrix (el, Trans, elmat);
220  elmat.Lump();
221 }
222 
223 void InverseIntegrator::SetIntRule(const IntegrationRule *ir)
224 {
225  IntRule = ir;
226  integrator->SetIntRule(ir);
227 }
228 
229 void InverseIntegrator::AssembleElementMatrix(
231 {
232  integrator->AssembleElementMatrix(el, Trans, elmat);
233  elmat.Invert();
234 }
235 
236 void SumIntegrator::SetIntRule(const IntegrationRule *ir)
237 {
238  IntRule = ir;
239  for (int i = 0; i < integrators.Size(); i++)
240  {
241  integrators[i]->SetIntRule(ir);
242  }
243 }
244 
245 void SumIntegrator::AssembleElementMatrix(
247 {
248  MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
249 
250  integrators[0]->AssembleElementMatrix(el, Trans, elmat);
251  for (int i = 1; i < integrators.Size(); i++)
252  {
253  integrators[i]->AssembleElementMatrix(el, Trans, elem_mat);
254  elmat += elem_mat;
255  }
256 }
257 
258 void SumIntegrator::AssembleElementMatrix2(
259  const FiniteElement &el1, const FiniteElement &el2,
261 {
262  MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
263 
264  integrators[0]->AssembleElementMatrix2(el1, el2, Trans, elmat);
265  for (int i = 1; i < integrators.Size(); i++)
266  {
267  integrators[i]->AssembleElementMatrix2(el1, el2, Trans, elem_mat);
268  elmat += elem_mat;
269  }
270 }
271 
272 void SumIntegrator::AssembleFaceMatrix(
273  const FiniteElement &el1, const FiniteElement &el2,
275 {
276  MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
277 
278  integrators[0]->AssembleFaceMatrix(el1, el2, Trans, elmat);
279  for (int i = 1; i < integrators.Size(); i++)
280  {
281  integrators[i]->AssembleFaceMatrix(el1, el2, Trans, elem_mat);
282  elmat += elem_mat;
283  }
284 }
285 
286 void SumIntegrator::AssembleFaceMatrix(
287  const FiniteElement &tr_fe,
288  const FiniteElement &te_fe1, const FiniteElement &te_fe2,
290 {
291  MFEM_ASSERT(integrators.Size() > 0, "empty SumIntegrator.");
292 
293  integrators[0]->AssembleFaceMatrix(tr_fe, te_fe1, te_fe2, Trans, elmat);
294  for (int i = 1; i < integrators.Size(); i++)
295  {
296  integrators[i]->AssembleFaceMatrix(tr_fe, te_fe1, te_fe2, Trans, elem_mat);
297  elmat += elem_mat;
298  }
299 }
300 
301 void SumIntegrator::AssemblePA(const FiniteElementSpace& fes)
302 {
303  for (int i = 0; i < integrators.Size(); i++)
304  {
305  integrators[i]->AssemblePA(fes);
306  }
307 }
308 
309 void SumIntegrator::AssembleDiagonalPA(Vector &diag)
310 {
311  for (int i = 0; i < integrators.Size(); i++)
312  {
313  integrators[i]->AssembleDiagonalPA(diag);
314  }
315 }
316 
317 void SumIntegrator::AssemblePAInteriorFaces(const FiniteElementSpace &fes)
318 {
319  for (int i = 0; i < integrators.Size(); i++)
320  {
321  integrators[i]->AssemblePAInteriorFaces(fes);
322  }
323 }
324 
325 void SumIntegrator::AssemblePABoundaryFaces(const FiniteElementSpace &fes)
326 {
327  for (int i = 0; i < integrators.Size(); i++)
328  {
329  integrators[i]->AssemblePABoundaryFaces(fes);
330  }
331 }
332 
333 void SumIntegrator::AddMultPA(const Vector& x, Vector& y) const
334 {
335  for (int i = 0; i < integrators.Size(); i++)
336  {
337  integrators[i]->AddMultPA(x, y);
338  }
339 }
340 
341 void SumIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
342 {
343  for (int i = 0; i < integrators.Size(); i++)
344  {
345  integrators[i]->AddMultTransposePA(x, y);
346  }
347 }
348 
349 void SumIntegrator::AssembleMF(const FiniteElementSpace &fes)
350 {
351  for (int i = 0; i < integrators.Size(); i++)
352  {
353  integrators[i]->AssembleMF(fes);
354  }
355 }
356 
357 void SumIntegrator::AddMultMF(const Vector& x, Vector& y) const
358 {
359  for (int i = 0; i < integrators.Size(); i++)
360  {
361  integrators[i]->AddMultTransposeMF(x, y);
362  }
363 }
364 
365 void SumIntegrator::AddMultTransposeMF(const Vector &x, Vector &y) const
366 {
367  for (int i = 0; i < integrators.Size(); i++)
368  {
369  integrators[i]->AddMultMF(x, y);
370  }
371 }
372 
373 void SumIntegrator::AssembleDiagonalMF(Vector &diag)
374 {
375  for (int i = 0; i < integrators.Size(); i++)
376  {
377  integrators[i]->AssembleDiagonalMF(diag);
378  }
379 }
380 
381 void SumIntegrator::AssembleEA(const FiniteElementSpace &fes, Vector &emat,
382  const bool add)
383 {
384  for (int i = 0; i < integrators.Size(); i++)
385  {
386  integrators[i]->AssembleEA(fes, emat, add);
387  }
388 }
389 
390 void SumIntegrator::AssembleEAInteriorFaces(const FiniteElementSpace &fes,
391  Vector &ea_data_int,
392  Vector &ea_data_ext,
393  const bool add)
394 {
395  for (int i = 0; i < integrators.Size(); i++)
396  {
397  integrators[i]->AssembleEAInteriorFaces(fes,ea_data_int,ea_data_ext,add);
398  }
399 }
400 
401 void SumIntegrator::AssembleEABoundaryFaces(const FiniteElementSpace &fes,
402  Vector &ea_data_bdr,
403  const bool add)
404 {
405  for (int i = 0; i < integrators.Size(); i++)
406  {
407  integrators[i]->AssembleEABoundaryFaces(fes, ea_data_bdr, add);
408  }
409 }
410 
411 SumIntegrator::~SumIntegrator()
412 {
413  if (own_integrators)
414  {
415  for (int i = 0; i < integrators.Size(); i++)
416  {
417  delete integrators[i];
418  }
419  }
420 }
421 
422 void MixedScalarIntegrator::AssembleElementMatrix2(
423  const FiniteElement &trial_fe, const FiniteElement &test_fe,
425 {
426  MFEM_ASSERT(this->VerifyFiniteElementTypes(trial_fe, test_fe),
427  this->FiniteElementTypeFailureMessage());
428 
429  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
430  bool same_shapes = same_calc_shape && (&trial_fe == &test_fe);
431 
432 #ifdef MFEM_THREAD_SAFE
433  Vector test_shape(test_nd);
434  Vector trial_shape;
435 #else
436  test_shape.SetSize(test_nd);
437 #endif
438  if (same_shapes)
439  {
440  trial_shape.NewDataAndSize(test_shape.GetData(), trial_nd);
441  }
442  else
443  {
444  trial_shape.SetSize(trial_nd);
445  }
446 
447  elmat.SetSize(test_nd, trial_nd);
448 
449  const IntegrationRule *ir = IntRule;
450  if (ir == NULL)
451  {
452  int ir_order = this->GetIntegrationOrder(trial_fe, test_fe, Trans);
453  ir = &IntRules.Get(trial_fe.GetGeomType(), ir_order);
454  }
455 
456  elmat = 0.0;
457  for (i = 0; i < ir->GetNPoints(); i++)
458  {
459  const IntegrationPoint &ip = ir->IntPoint(i);
460  Trans.SetIntPoint(&ip);
461 
462  this->CalcTestShape(test_fe, Trans, test_shape);
463  this->CalcTrialShape(trial_fe, Trans, trial_shape);
464 
465  double w = Trans.Weight() * ip.weight;
466 
467  if (Q)
468  {
469  w *= Q->Eval(Trans, ip);
470  }
471  AddMult_a_VWt(w, test_shape, trial_shape, elmat);
472  }
473 #ifndef MFEM_THREAD_SAFE
474  if (same_shapes)
475  {
476  trial_shape.SetDataAndSize(NULL, 0);
477  }
478 #endif
479 }
480 
481 void MixedVectorIntegrator::AssembleElementMatrix2(
482  const FiniteElement &trial_fe, const FiniteElement &test_fe,
484 {
485  MFEM_ASSERT(this->VerifyFiniteElementTypes(trial_fe, test_fe),
486  this->FiniteElementTypeFailureMessage());
487 
488  space_dim = Trans.GetSpaceDim();
489  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
490  int test_vdim = GetTestVDim(test_fe);
491  int trial_vdim = GetTrialVDim(trial_fe);
492  bool same_shapes = same_calc_shape && (&trial_fe == &test_fe);
493 
494  if (MQ)
495  {
496  MFEM_VERIFY(MQ->GetHeight() == test_vdim,
497  "Dimension mismatch in height of matrix coefficient.");
498  MFEM_VERIFY(MQ->GetWidth() == trial_vdim,
499  "Dimension mismatch in width of matrix coefficient.");
500  }
501  if (DQ)
502  {
503  MFEM_VERIFY(trial_vdim == test_vdim,
504  "Diagonal matrix coefficient requires matching "
505  "test and trial vector dimensions.");
506  MFEM_VERIFY(DQ->GetVDim() == trial_vdim,
507  "Dimension mismatch in diagonal matrix coefficient.");
508  }
509  if (VQ)
510  {
511  MFEM_VERIFY(VQ->GetVDim() == 3, "Vector coefficient must have "
512  "dimension equal to three.");
513  }
514 
515 #ifdef MFEM_THREAD_SAFE
516  Vector V(VQ ? VQ->GetVDim() : 0);
517  Vector D(DQ ? DQ->GetVDim() : 0);
518  DenseMatrix M(MQ ? MQ->GetHeight() : 0, MQ ? MQ->GetWidth() : 0);
519  DenseMatrix test_shape(test_nd, test_vdim);
520  DenseMatrix trial_shape;
521  DenseMatrix shape_tmp(test_nd, trial_vdim);
522 #else
523  V.SetSize(VQ ? VQ->GetVDim() : 0);
524  D.SetSize(DQ ? DQ->GetVDim() : 0);
525  M.SetSize(MQ ? MQ->GetHeight() : 0, MQ ? MQ->GetWidth() : 0);
526  test_shape.SetSize(test_nd, test_vdim);
527  shape_tmp.SetSize(test_nd, trial_vdim);
528 #endif
529  if (same_shapes)
530  {
531  trial_shape.Reset(test_shape.Data(), trial_nd, trial_vdim);
532  }
533  else
534  {
535  trial_shape.SetSize(trial_nd, trial_vdim);
536  }
537 
538  elmat.SetSize(test_nd, trial_nd);
539 
540  const IntegrationRule *ir = IntRule;
541  if (ir == NULL)
542  {
543  int ir_order = this->GetIntegrationOrder(trial_fe, test_fe, Trans);
544  ir = &IntRules.Get(trial_fe.GetGeomType(), ir_order);
545  }
546 
547  elmat = 0.0;
548  for (i = 0; i < ir->GetNPoints(); i++)
549  {
550  const IntegrationPoint &ip = ir->IntPoint(i);
551  Trans.SetIntPoint(&ip);
552 
553  this->CalcTestShape(test_fe, Trans, test_shape);
554  if (!same_shapes)
555  {
556  this->CalcTrialShape(trial_fe, Trans, trial_shape);
557  }
558 
559  double w = Trans.Weight() * ip.weight;
560 
561  if (MQ)
562  {
563  MQ->Eval(M, Trans, ip);
564  M *= w;
565  Mult(test_shape, M, shape_tmp);
566  AddMultABt(shape_tmp, trial_shape, elmat);
567  }
568  else if (DQ)
569  {
570  DQ->Eval(D, Trans, ip);
571  D *= w;
572  AddMultADBt(test_shape, D, trial_shape, elmat);
573  }
574  else if (VQ)
575  {
576  VQ->Eval(V, Trans, ip);
577  V *= w;
578 
579  for (int j=0; j<test_nd; j++)
580  {
581  // Compute shape_tmp = test_shape x V
582  // V will always be of length 3
583  // shape_dim and test_shape could have reduced dimension
584  // i.e. 1D or 2D
585  if (test_vdim == 3 && trial_vdim == 3)
586  {
587  shape_tmp(j,0) = test_shape(j,1) * V(2) -
588  test_shape(j,2) * V(1);
589  shape_tmp(j,1) = test_shape(j,2) * V(0) -
590  test_shape(j,0) * V(2);
591  shape_tmp(j,2) = test_shape(j,0) * V(1) -
592  test_shape(j,1) * V(0);
593  }
594  else if (test_vdim == 3 && trial_vdim == 2)
595  {
596  shape_tmp(j,0) = test_shape(j,1) * V(2) -
597  test_shape(j,2) * V(1);
598  shape_tmp(j,1) = test_shape(j,2) * V(0) -
599  test_shape(j,0) * V(2);
600  }
601  else if (test_vdim == 3 && trial_vdim == 1)
602  {
603  shape_tmp(j,0) = test_shape(j,1) * V(2) -
604  test_shape(j,2) * V(1);
605  }
606  else if (test_vdim == 2 && trial_vdim == 3)
607  {
608  shape_tmp(j,0) = test_shape(j,1) * V(2);
609  shape_tmp(j,1) = -test_shape(j,0) * V(2);
610  shape_tmp(j,2) = test_shape(j,0) * V(1) -
611  test_shape(j,1) * V(0);
612  }
613  else if (test_vdim == 2 && trial_vdim == 2)
614  {
615  shape_tmp(j,0) = test_shape(j,1) * V(2);
616  shape_tmp(j,1) = -test_shape(j,0) * V(2);
617  }
618  else if (test_vdim == 1 && trial_vdim == 3)
619  {
620  shape_tmp(j,0) = 0.0;
621  shape_tmp(j,1) = -test_shape(j,0) * V(2);
622  shape_tmp(j,2) = test_shape(j,0) * V(1);
623  }
624  else if (test_vdim == 1 && trial_vdim == 1)
625  {
626  shape_tmp(j,0) = 0.0;
627  }
628  }
629  AddMultABt(shape_tmp, trial_shape, elmat);
630  }
631  else
632  {
633  if (Q)
634  {
635  w *= Q -> Eval (Trans, ip);
636  }
637  if (same_shapes)
638  {
639  AddMult_a_AAt (w, test_shape, elmat);
640  }
641  else
642  {
643  AddMult_a_ABt (w, test_shape, trial_shape, elmat);
644  }
645  }
646  }
647 #ifndef MFEM_THREAD_SAFE
648  if (same_shapes)
649  {
650  trial_shape.ClearExternalData();
651  }
652 #endif
653 }
654 
655 void MixedScalarVectorIntegrator::AssembleElementMatrix2(
656  const FiniteElement &trial_fe, const FiniteElement &test_fe,
658 {
659  MFEM_ASSERT(this->VerifyFiniteElementTypes(trial_fe, test_fe),
660  this->FiniteElementTypeFailureMessage());
661 
662  MFEM_VERIFY(VQ, "MixedScalarVectorIntegrator: "
663  "VectorCoefficient must be set");
664 
665  const FiniteElement * vec_fe = transpose?&trial_fe:&test_fe;
666  const FiniteElement * sca_fe = transpose?&test_fe:&trial_fe;
667 
668  space_dim = Trans.GetSpaceDim();
669  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
670  int sca_nd = sca_fe->GetDof();
671  int vec_nd = vec_fe->GetDof();
672  int vdim = GetVDim(*vec_fe);
673  double vtmp;
674 
675  MFEM_VERIFY(VQ->GetVDim() == vdim, "MixedScalarVectorIntegrator: "
676  "Dimensions of VectorCoefficient and Vector-valued basis "
677  "functions must match");
678 
679 #ifdef MFEM_THREAD_SAFE
680  Vector V(vdim);
681  DenseMatrix vshape(vec_nd, vdim);
682  Vector shape(sca_nd);
683  Vector vshape_tmp(vec_nd);
684 #else
685  V.SetSize(vdim);
686  vshape.SetSize(vec_nd, vdim);
687  shape.SetSize(sca_nd);
688  vshape_tmp.SetSize(vec_nd);
689 #endif
690 
691  Vector V_test(transpose?shape.GetData():vshape_tmp.GetData(),test_nd);
692  Vector W_trial(transpose?vshape_tmp.GetData():shape.GetData(),trial_nd);
693 
694  elmat.SetSize(test_nd, trial_nd);
695 
696  const IntegrationRule *ir = IntRule;
697  if (ir == NULL)
698  {
699  int ir_order = this->GetIntegrationOrder(trial_fe, test_fe, Trans);
700  ir = &IntRules.Get(trial_fe.GetGeomType(), ir_order);
701  }
702 
703  elmat = 0.0;
704  for (i = 0; i < ir->GetNPoints(); i++)
705  {
706  const IntegrationPoint &ip = ir->IntPoint(i);
707  Trans.SetIntPoint(&ip);
708 
709  this->CalcShape(*sca_fe, Trans, shape);
710  this->CalcVShape(*vec_fe, Trans, vshape);
711 
712  double w = Trans.Weight() * ip.weight;
713 
714  VQ->Eval(V, Trans, ip);
715  V *= w;
716 
717  if ( vdim == 2 && cross_2d )
718  {
719  vtmp = V[0];
720  V[0] = -V[1];
721  V[1] = vtmp;
722  }
723 
724  vshape.Mult(V,vshape_tmp);
725  AddMultVWt(V_test, W_trial, elmat);
726  }
727 }
728 
729 
730 void GradientIntegrator::AssembleElementMatrix2(
731  const FiniteElement &trial_fe, const FiniteElement &test_fe,
733 {
734  dim = test_fe.GetDim();
735  int trial_dof = trial_fe.GetDof();
736  int test_dof = test_fe.GetDof();
737  double c;
738  Vector d_col;
739 
740  dshape.SetSize(trial_dof, dim);
741  gshape.SetSize(trial_dof, dim);
742  Jadj.SetSize(dim);
743  shape.SetSize(test_dof);
744  elmat.SetSize(dim * test_dof, trial_dof);
745 
746  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe,
747  Trans);
748 
749  elmat = 0.0;
750  elmat_comp.SetSize(test_dof, trial_dof);
751 
752  for (int i = 0; i < ir->GetNPoints(); i++)
753  {
754  const IntegrationPoint &ip = ir->IntPoint(i);
755 
756  trial_fe.CalcDShape(ip, dshape);
757  test_fe.CalcShape(ip, shape);
758 
759  Trans.SetIntPoint(&ip);
760  CalcAdjugate(Trans.Jacobian(), Jadj);
761 
762  Mult(dshape, Jadj, gshape);
763 
764  c = ip.weight;
765  if (Q)
766  {
767  c *= Q->Eval(Trans, ip);
768  }
769  shape *= c;
770 
771  for (int d = 0; d < dim; ++d)
772  {
773  gshape.GetColumnReference(d, d_col);
774  MultVWt(shape, d_col, elmat_comp);
775  for (int jj = 0; jj < trial_dof; ++jj)
776  {
777  for (int ii = 0; ii < test_dof; ++ii)
778  {
779  elmat(d * test_dof + ii, jj) += elmat_comp(ii, jj);
780  }
781  }
782  }
783  }
784 }
785 
787  &trial_fe,
788  const FiniteElement &test_fe,
790 {
791  int order = Trans.OrderGrad(&trial_fe) + test_fe.GetOrder() + Trans.OrderJ();
792  return IntRules.Get(trial_fe.GetGeomType(), order);
793 }
794 
795 
796 void DiffusionIntegrator::AssembleElementMatrix
798  DenseMatrix &elmat )
799 {
800  int nd = el.GetDof();
801  dim = el.GetDim();
802  int spaceDim = Trans.GetSpaceDim();
803  bool square = (dim == spaceDim);
804  double w;
805 
806  if (VQ)
807  {
808  MFEM_VERIFY(VQ->GetVDim() == spaceDim,
809  "Unexpected dimension for VectorCoefficient");
810  }
811  if (MQ)
812  {
813  MFEM_VERIFY(MQ->GetWidth() == spaceDim,
814  "Unexpected width for MatrixCoefficient");
815  MFEM_VERIFY(MQ->GetHeight() == spaceDim,
816  "Unexpected height for MatrixCoefficient");
817  }
818 
819 #ifdef MFEM_THREAD_SAFE
820  DenseMatrix dshape(nd, dim), dshapedxt(nd, spaceDim);
821  DenseMatrix dshapedxt_m(nd, MQ ? spaceDim : 0);
822  DenseMatrix M(MQ ? spaceDim : 0);
823  Vector D(VQ ? VQ->GetVDim() : 0);
824 #else
825  dshape.SetSize(nd, dim);
826  dshapedxt.SetSize(nd, spaceDim);
827  dshapedxt_m.SetSize(nd, MQ ? spaceDim : 0);
828  M.SetSize(MQ ? spaceDim : 0);
829  D.SetSize(VQ ? VQ->GetVDim() : 0);
830 #endif
831  elmat.SetSize(nd);
832 
833  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el);
834 
835  elmat = 0.0;
836  for (int i = 0; i < ir->GetNPoints(); i++)
837  {
838  const IntegrationPoint &ip = ir->IntPoint(i);
839  el.CalcDShape(ip, dshape);
840 
841  Trans.SetIntPoint(&ip);
842  w = Trans.Weight();
843  w = ip.weight / (square ? w : w*w*w);
844  // AdjugateJacobian = / adj(J), if J is square
845  // \ adj(J^t.J).J^t, otherwise
846  Mult(dshape, Trans.AdjugateJacobian(), dshapedxt);
847  if (MQ)
848  {
849  MQ->Eval(M, Trans, ip);
850  M *= w;
851  Mult(dshapedxt, M, dshapedxt_m);
852  AddMultABt(dshapedxt_m, dshapedxt, elmat);
853  }
854  else if (VQ)
855  {
856  VQ->Eval(D, Trans, ip);
857  D *= w;
858  AddMultADAt(dshapedxt, D, elmat);
859  }
860  else
861  {
862  if (Q)
863  {
864  w *= Q->Eval(Trans, ip);
865  }
866  AddMult_a_AAt(w, dshapedxt, elmat);
867  }
868  }
869 }
870 
871 void DiffusionIntegrator::AssembleElementMatrix2(
872  const FiniteElement &trial_fe, const FiniteElement &test_fe,
874 {
875  int tr_nd = trial_fe.GetDof();
876  int te_nd = test_fe.GetDof();
877  dim = trial_fe.GetDim();
878  int spaceDim = Trans.GetSpaceDim();
879  bool square = (dim == spaceDim);
880  double w;
881 
882  if (VQ)
883  {
884  MFEM_VERIFY(VQ->GetVDim() == spaceDim,
885  "Unexpected dimension for VectorCoefficient");
886  }
887  if (MQ)
888  {
889  MFEM_VERIFY(MQ->GetWidth() == spaceDim,
890  "Unexpected width for MatrixCoefficient");
891  MFEM_VERIFY(MQ->GetHeight() == spaceDim,
892  "Unexpected height for MatrixCoefficient");
893  }
894 
895 #ifdef MFEM_THREAD_SAFE
896  DenseMatrix dshape(tr_nd, dim), dshapedxt(tr_nd, spaceDim);
897  DenseMatrix te_dshape(te_nd, dim), te_dshapedxt(te_nd, spaceDim);
898  DenseMatrix invdfdx(dim, spaceDim);
899  DenseMatrix dshapedxt_m(te_nd, MQ ? spaceDim : 0);
900  DenseMatrix M(MQ ? spaceDim : 0);
901  Vector D(VQ ? VQ->GetVDim() : 0);
902 #else
903  dshape.SetSize(tr_nd, dim);
904  dshapedxt.SetSize(tr_nd, spaceDim);
905  te_dshape.SetSize(te_nd, dim);
906  te_dshapedxt.SetSize(te_nd, spaceDim);
907  invdfdx.SetSize(dim, spaceDim);
908  dshapedxt_m.SetSize(te_nd, MQ ? spaceDim : 0);
909  M.SetSize(MQ ? spaceDim : 0);
910  D.SetSize(VQ ? VQ->GetVDim() : 0);
911 #endif
912  elmat.SetSize(te_nd, tr_nd);
913 
914  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe);
915 
916  elmat = 0.0;
917  for (int i = 0; i < ir->GetNPoints(); i++)
918  {
919  const IntegrationPoint &ip = ir->IntPoint(i);
920  trial_fe.CalcDShape(ip, dshape);
921  test_fe.CalcDShape(ip, te_dshape);
922 
923  Trans.SetIntPoint(&ip);
924  CalcAdjugate(Trans.Jacobian(), invdfdx);
925  w = Trans.Weight();
926  w = ip.weight / (square ? w : w*w*w);
927  Mult(dshape, invdfdx, dshapedxt);
928  Mult(te_dshape, invdfdx, te_dshapedxt);
929  // invdfdx, dshape, and te_dshape no longer needed
930  if (MQ)
931  {
932  MQ->Eval(M, Trans, ip);
933  M *= w;
934  Mult(te_dshapedxt, M, dshapedxt_m);
935  AddMultABt(dshapedxt_m, dshapedxt, elmat);
936  }
937  else if (VQ)
938  {
939  VQ->Eval(D, Trans, ip);
940  D *= w;
941  AddMultADAt(dshapedxt, D, elmat);
942  }
943  else
944  {
945  if (Q)
946  {
947  w *= Q->Eval(Trans, ip);
948  }
949  dshapedxt *= w;
950  AddMultABt(te_dshapedxt, dshapedxt, elmat);
951  }
952  }
953 }
954 
955 void DiffusionIntegrator::AssembleElementVector(
956  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
957  Vector &elvect)
958 {
959  int nd = el.GetDof();
960  dim = el.GetDim();
961  int spaceDim = Tr.GetSpaceDim();
962  double w;
963 
964  if (VQ)
965  {
966  MFEM_VERIFY(VQ->GetVDim() == spaceDim,
967  "Unexpected dimension for VectorCoefficient");
968  }
969  if (MQ)
970  {
971  MFEM_VERIFY(MQ->GetWidth() == spaceDim,
972  "Unexpected width for MatrixCoefficient");
973  MFEM_VERIFY(MQ->GetHeight() == spaceDim,
974  "Unexpected height for MatrixCoefficient");
975  }
976 
977 #ifdef MFEM_THREAD_SAFE
978  DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim), M(MQ ? spaceDim : 0);
979  Vector D(VQ ? VQ->GetVDim() : 0);
980 #else
981  dshape.SetSize(nd,dim);
982  invdfdx.SetSize(dim, spaceDim);
983  M.SetSize(MQ ? spaceDim : 0);
984  D.SetSize(VQ ? VQ->GetVDim() : 0);
985 #endif
986  vec.SetSize(dim);
987  vecdxt.SetSize((VQ || MQ) ? spaceDim : 0);
988  pointflux.SetSize(spaceDim);
989 
990  elvect.SetSize(nd);
991 
992  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el);
993 
994  elvect = 0.0;
995  for (int i = 0; i < ir->GetNPoints(); i++)
996  {
997  const IntegrationPoint &ip = ir->IntPoint(i);
998  el.CalcDShape(ip, dshape);
999 
1000  Tr.SetIntPoint(&ip);
1001  CalcAdjugate(Tr.Jacobian(), invdfdx); // invdfdx = adj(J)
1002  w = ip.weight / Tr.Weight();
1003 
1004  if (!MQ && !VQ)
1005  {
1006  dshape.MultTranspose(elfun, vec);
1007  invdfdx.MultTranspose(vec, pointflux);
1008  if (Q)
1009  {
1010  w *= Q->Eval(Tr, ip);
1011  }
1012  }
1013  else
1014  {
1015  dshape.MultTranspose(elfun, vec);
1016  invdfdx.MultTranspose(vec, vecdxt);
1017  if (MQ)
1018  {
1019  MQ->Eval(M, Tr, ip);
1020  M.Mult(vecdxt, pointflux);
1021  }
1022  else
1023  {
1024  VQ->Eval(D, Tr, ip);
1025  for (int j=0; j<spaceDim; ++j)
1026  {
1027  pointflux[j] = D[j] * vecdxt[j];
1028  }
1029  }
1030  }
1031  pointflux *= w;
1032  invdfdx.Mult(pointflux, vec);
1033  dshape.AddMult(vec, elvect);
1034  }
1035 }
1036 
1037 void DiffusionIntegrator::ComputeElementFlux
1039  Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef,
1040  const IntegrationRule *ir)
1041 {
1042  int nd, spaceDim, fnd;
1043 
1044  nd = el.GetDof();
1045  dim = el.GetDim();
1046  spaceDim = Trans.GetSpaceDim();
1047 
1048  if (VQ)
1049  {
1050  MFEM_VERIFY(VQ->GetVDim() == spaceDim,
1051  "Unexpected dimension for VectorCoefficient");
1052  }
1053  if (MQ)
1054  {
1055  MFEM_VERIFY(MQ->GetWidth() == spaceDim,
1056  "Unexpected width for MatrixCoefficient");
1057  MFEM_VERIFY(MQ->GetHeight() == spaceDim,
1058  "Unexpected height for MatrixCoefficient");
1059  }
1060 
1061 #ifdef MFEM_THREAD_SAFE
1062  DenseMatrix dshape(nd,dim), invdfdx(dim, spaceDim);
1063  DenseMatrix M(MQ ? spaceDim : 0);
1064  Vector D(VQ ? VQ->GetVDim() : 0);
1065 #else
1066  dshape.SetSize(nd,dim);
1067  invdfdx.SetSize(dim, spaceDim);
1068  M.SetSize(MQ ? spaceDim : 0);
1069  D.SetSize(VQ ? VQ->GetVDim() : 0);
1070 #endif
1071  vec.SetSize(dim);
1072  vecdxt.SetSize(spaceDim);
1073  pointflux.SetSize(MQ || VQ ? spaceDim : 0);
1074 
1075  if (!ir)
1076  {
1077  ir = &fluxelem.GetNodes();
1078  }
1079  fnd = ir->GetNPoints();
1080  flux.SetSize( fnd * spaceDim );
1081 
1082  for (int i = 0; i < fnd; i++)
1083  {
1084  const IntegrationPoint &ip = ir->IntPoint(i);
1085  el.CalcDShape(ip, dshape);
1086  dshape.MultTranspose(u, vec);
1087 
1088  Trans.SetIntPoint (&ip);
1089  CalcInverse(Trans.Jacobian(), invdfdx);
1090  invdfdx.MultTranspose(vec, vecdxt);
1091 
1092  if (with_coef)
1093  {
1094  if (!MQ && !VQ)
1095  {
1096  if (Q)
1097  {
1098  vecdxt *= Q->Eval(Trans,ip);
1099  }
1100  for (int j = 0; j < spaceDim; j++)
1101  {
1102  flux(fnd*j+i) = vecdxt(j);
1103  }
1104  }
1105  else
1106  {
1107  if (MQ)
1108  {
1109  MQ->Eval(M, Trans, ip);
1110  M.Mult(vecdxt, pointflux);
1111  }
1112  else
1113  {
1114  VQ->Eval(D, Trans, ip);
1115  for (int j=0; j<spaceDim; ++j)
1116  {
1117  pointflux[j] = D[j] * vecdxt[j];
1118  }
1119  }
1120  for (int j = 0; j < spaceDim; j++)
1121  {
1122  flux(fnd*j+i) = pointflux(j);
1123  }
1124  }
1125  }
1126  else
1127  {
1128  for (int j = 0; j < spaceDim; j++)
1129  {
1130  flux(fnd*j+i) = vecdxt(j);
1131  }
1132  }
1133  }
1134 }
1135 
1136 double DiffusionIntegrator::ComputeFluxEnergy
1138  Vector &flux, Vector* d_energy)
1139 {
1140  int nd = fluxelem.GetDof();
1141  dim = fluxelem.GetDim();
1142  int spaceDim = Trans.GetSpaceDim();
1143 
1144 #ifdef MFEM_THREAD_SAFE
1145  DenseMatrix M;
1146  Vector D(VQ ? VQ->GetVDim() : 0);
1147 #else
1148  D.SetSize(VQ ? VQ->GetVDim() : 0);
1149 #endif
1150 
1151  shape.SetSize(nd);
1152  pointflux.SetSize(spaceDim);
1153  if (d_energy) { vec.SetSize(spaceDim); }
1154  if (MQ) { M.SetSize(spaceDim); }
1155 
1156  int order = 2 * fluxelem.GetOrder(); // <--
1157  const IntegrationRule *ir = &IntRules.Get(fluxelem.GetGeomType(), order);
1158 
1159  double energy = 0.0;
1160  if (d_energy) { *d_energy = 0.0; }
1161 
1162  for (int i = 0; i < ir->GetNPoints(); i++)
1163  {
1164  const IntegrationPoint &ip = ir->IntPoint(i);
1165  fluxelem.CalcShape(ip, shape);
1166 
1167  pointflux = 0.0;
1168  for (int k = 0; k < spaceDim; k++)
1169  {
1170  for (int j = 0; j < nd; j++)
1171  {
1172  pointflux(k) += flux(k*nd+j)*shape(j);
1173  }
1174  }
1175 
1176  Trans.SetIntPoint(&ip);
1177  double w = Trans.Weight() * ip.weight;
1178 
1179  if (MQ)
1180  {
1181  MQ->Eval(M, Trans, ip);
1182  energy += w * M.InnerProduct(pointflux, pointflux);
1183  }
1184  else if (VQ)
1185  {
1186  VQ->Eval(D, Trans, ip);
1187  D *= pointflux;
1188  energy += w * (D * pointflux);
1189  }
1190  else
1191  {
1192  double e = (pointflux * pointflux);
1193  if (Q) { e *= Q->Eval(Trans, ip); }
1194  energy += w * e;
1195  }
1196 
1197  if (d_energy)
1198  {
1199  // transform pointflux to the ref. domain and integrate the components
1200  Trans.Jacobian().MultTranspose(pointflux, vec);
1201  for (int k = 0; k < dim; k++)
1202  {
1203  (*d_energy)[k] += w * vec[k] * vec[k];
1204  }
1205  // TODO: Q, VQ, MQ
1206  }
1207  }
1208 
1209  return energy;
1210 }
1211 
1213  const FiniteElement &trial_fe, const FiniteElement &test_fe)
1214 {
1215  int order;
1216  if (trial_fe.Space() == FunctionSpace::Pk)
1217  {
1218  order = trial_fe.GetOrder() + test_fe.GetOrder() - 2;
1219  }
1220  else
1221  {
1222  // order = 2*el.GetOrder() - 2; // <-- this seems to work fine too
1223  order = trial_fe.GetOrder() + test_fe.GetOrder() + trial_fe.GetDim() - 1;
1224  }
1225 
1226  if (trial_fe.Space() == FunctionSpace::rQk)
1227  {
1228  return RefinedIntRules.Get(trial_fe.GetGeomType(), order);
1229  }
1230  return IntRules.Get(trial_fe.GetGeomType(), order);
1231 }
1232 
1233 
1234 void MassIntegrator::AssembleElementMatrix
1236  DenseMatrix &elmat )
1237 {
1238  int nd = el.GetDof();
1239  // int dim = el.GetDim();
1240  double w;
1241 
1242 #ifdef MFEM_THREAD_SAFE
1243  Vector shape;
1244 #endif
1245  elmat.SetSize(nd);
1246  shape.SetSize(nd);
1247 
1248  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, Trans);
1249 
1250  elmat = 0.0;
1251  for (int i = 0; i < ir->GetNPoints(); i++)
1252  {
1253  const IntegrationPoint &ip = ir->IntPoint(i);
1254  Trans.SetIntPoint (&ip);
1255 
1256  el.CalcPhysShape(Trans, shape);
1257 
1258  w = Trans.Weight() * ip.weight;
1259  if (Q)
1260  {
1261  w *= Q -> Eval(Trans, ip);
1262  }
1263 
1264  AddMult_a_VVt(w, shape, elmat);
1265  }
1266 }
1267 
1268 void MassIntegrator::AssembleElementMatrix2(
1269  const FiniteElement &trial_fe, const FiniteElement &test_fe,
1271 {
1272  int tr_nd = trial_fe.GetDof();
1273  int te_nd = test_fe.GetDof();
1274  double w;
1275 
1276 #ifdef MFEM_THREAD_SAFE
1277  Vector shape, te_shape;
1278 #endif
1279  elmat.SetSize(te_nd, tr_nd);
1280  shape.SetSize(tr_nd);
1281  te_shape.SetSize(te_nd);
1282 
1283  const IntegrationRule *ir = IntRule ? IntRule :
1284  &GetRule(trial_fe, test_fe, Trans);
1285 
1286  elmat = 0.0;
1287  for (int i = 0; i < ir->GetNPoints(); i++)
1288  {
1289  const IntegrationPoint &ip = ir->IntPoint(i);
1290  trial_fe.CalcShape(ip, shape);
1291  test_fe.CalcShape(ip, te_shape);
1292 
1293  Trans.SetIntPoint (&ip);
1294  w = Trans.Weight() * ip.weight;
1295  if (Q)
1296  {
1297  w *= Q -> Eval(Trans, ip);
1298  }
1299 
1300  te_shape *= w;
1301  AddMultVWt(te_shape, shape, elmat);
1302  }
1303 }
1304 
1306  const FiniteElement &test_fe,
1308 {
1309  // int order = trial_fe.GetOrder() + test_fe.GetOrder();
1310  const int order = trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW();
1311 
1312  if (trial_fe.Space() == FunctionSpace::rQk)
1313  {
1314  return RefinedIntRules.Get(trial_fe.GetGeomType(), order);
1315  }
1316  return IntRules.Get(trial_fe.GetGeomType(), order);
1317 }
1318 
1319 
1320 void BoundaryMassIntegrator::AssembleFaceMatrix(
1321  const FiniteElement &el1, const FiniteElement &el2,
1323 {
1324  MFEM_ASSERT(Trans.Elem2No < 0,
1325  "support for interior faces is not implemented");
1326 
1327  int nd1 = el1.GetDof();
1328  double w;
1329 
1330 #ifdef MFEM_THREAD_SAFE
1331  Vector shape;
1332 #endif
1333  elmat.SetSize(nd1);
1334  shape.SetSize(nd1);
1335 
1336  const IntegrationRule *ir = IntRule;
1337  if (ir == NULL)
1338  {
1339  int order = 2 * el1.GetOrder();
1340 
1341  ir = &IntRules.Get(Trans.GetGeometryType(), order);
1342  }
1343 
1344  elmat = 0.0;
1345  for (int i = 0; i < ir->GetNPoints(); i++)
1346  {
1347  const IntegrationPoint &ip = ir->IntPoint(i);
1348 
1349  // Set the integration point in the face and the neighboring element
1350  Trans.SetAllIntPoints(&ip);
1351 
1352  // Access the neighboring element's integration point
1353  const IntegrationPoint &eip = Trans.GetElement1IntPoint();
1354  el1.CalcShape(eip, shape);
1355 
1356  w = Trans.Weight() * ip.weight;
1357  if (Q)
1358  {
1359  w *= Q -> Eval(Trans, ip);
1360  }
1361 
1362  AddMult_a_VVt(w, shape, elmat);
1363  }
1364 }
1365 
1366 void ConvectionIntegrator::AssembleElementMatrix(
1368 {
1369  int nd = el.GetDof();
1370  dim = el.GetDim();
1371 
1372 #ifdef MFEM_THREAD_SAFE
1373  DenseMatrix dshape, adjJ, Q_ir;
1374  Vector shape, vec2, BdFidxT;
1375 #endif
1376  elmat.SetSize(nd);
1377  dshape.SetSize(nd,dim);
1378  adjJ.SetSize(dim);
1379  shape.SetSize(nd);
1380  vec2.SetSize(dim);
1381  BdFidxT.SetSize(nd);
1382 
1383  Vector vec1;
1384 
1385  const IntegrationRule *ir = IntRule;
1386  if (ir == NULL)
1387  {
1388  int order = Trans.OrderGrad(&el) + Trans.Order() + el.GetOrder();
1389  ir = &IntRules.Get(el.GetGeomType(), order);
1390  }
1391 
1392  Q->Eval(Q_ir, Trans, *ir);
1393 
1394  elmat = 0.0;
1395  for (int i = 0; i < ir->GetNPoints(); i++)
1396  {
1397  const IntegrationPoint &ip = ir->IntPoint(i);
1398  el.CalcDShape(ip, dshape);
1399  el.CalcShape(ip, shape);
1400 
1401  Trans.SetIntPoint(&ip);
1402  CalcAdjugate(Trans.Jacobian(), adjJ);
1403  Q_ir.GetColumnReference(i, vec1);
1404  vec1 *= alpha * ip.weight;
1405 
1406  adjJ.Mult(vec1, vec2);
1407  dshape.Mult(vec2, BdFidxT);
1408 
1409  AddMultVWt(shape, BdFidxT, elmat);
1410  }
1411 }
1412 
1413 
1414 void GroupConvectionIntegrator::AssembleElementMatrix(
1416 {
1417  int nd = el.GetDof();
1418  int dim = el.GetDim();
1419 
1420  elmat.SetSize(nd);
1421  dshape.SetSize(nd,dim);
1422  adjJ.SetSize(dim);
1423  shape.SetSize(nd);
1424  grad.SetSize(nd,dim);
1425 
1426  const IntegrationRule *ir = IntRule;
1427  if (ir == NULL)
1428  {
1429  int order = Trans.OrderGrad(&el) + el.GetOrder();
1430  ir = &IntRules.Get(el.GetGeomType(), order);
1431  }
1432 
1433  Q->Eval(Q_nodal, Trans, el.GetNodes()); // sets the size of Q_nodal
1434 
1435  elmat = 0.0;
1436  for (int i = 0; i < ir->GetNPoints(); i++)
1437  {
1438  const IntegrationPoint &ip = ir->IntPoint(i);
1439  el.CalcDShape(ip, dshape);
1440  el.CalcShape(ip, shape);
1441 
1442  Trans.SetIntPoint(&ip);
1443  CalcAdjugate(Trans.Jacobian(), adjJ);
1444 
1445  Mult(dshape, adjJ, grad);
1446 
1447  double w = alpha * ip.weight;
1448 
1449  // elmat(k,l) += \sum_s w*shape(k)*Q_nodal(s,k)*grad(l,s)
1450  for (int k = 0; k < nd; k++)
1451  {
1452  double wsk = w*shape(k);
1453  for (int l = 0; l < nd; l++)
1454  {
1455  double a = 0.0;
1456  for (int s = 0; s < dim; s++)
1457  {
1458  a += Q_nodal(s,k)*grad(l,s);
1459  }
1460  elmat(k,l) += wsk*a;
1461  }
1462  }
1463  }
1464 }
1465 
1467  const FiniteElement &trial_fe, const FiniteElement &test_fe,
1469 {
1470  int order = Trans.OrderGrad(&trial_fe) + Trans.Order() + test_fe.GetOrder();
1471 
1472  return IntRules.Get(trial_fe.GetGeomType(), order);
1473 }
1474 
1477 {
1478  return GetRule(el,el,Trans);
1479 }
1480 
1481 void VectorMassIntegrator::AssembleElementMatrix
1483  DenseMatrix &elmat )
1484 {
1485  int nd = el.GetDof();
1486  int spaceDim = Trans.GetSpaceDim();
1487 
1488  double norm;
1489 
1490  // If vdim is not set, set it to the space dimension
1491  vdim = (vdim == -1) ? spaceDim : vdim;
1492 
1493  elmat.SetSize(nd*vdim);
1494  shape.SetSize(nd);
1495  partelmat.SetSize(nd);
1496  if (VQ)
1497  {
1498  vec.SetSize(vdim);
1499  }
1500  else if (MQ)
1501  {
1502  mcoeff.SetSize(vdim);
1503  }
1504 
1505  const IntegrationRule *ir = IntRule;
1506  if (ir == NULL)
1507  {
1508  int order = 2 * el.GetOrder() + Trans.OrderW() + Q_order;
1509 
1510  if (el.Space() == FunctionSpace::rQk)
1511  {
1512  ir = &RefinedIntRules.Get(el.GetGeomType(), order);
1513  }
1514  else
1515  {
1516  ir = &IntRules.Get(el.GetGeomType(), order);
1517  }
1518  }
1519 
1520  elmat = 0.0;
1521  for (int s = 0; s < ir->GetNPoints(); s++)
1522  {
1523  const IntegrationPoint &ip = ir->IntPoint(s);
1524  el.CalcShape(ip, shape);
1525 
1526  Trans.SetIntPoint (&ip);
1527  norm = ip.weight * Trans.Weight();
1528 
1529  MultVVt(shape, partelmat);
1530 
1531  if (VQ)
1532  {
1533  VQ->Eval(vec, Trans, ip);
1534  for (int k = 0; k < vdim; k++)
1535  {
1536  elmat.AddMatrix(norm*vec(k), partelmat, nd*k, nd*k);
1537  }
1538  }
1539  else if (MQ)
1540  {
1541  MQ->Eval(mcoeff, Trans, ip);
1542  for (int i = 0; i < vdim; i++)
1543  for (int j = 0; j < vdim; j++)
1544  {
1545  elmat.AddMatrix(norm*mcoeff(i,j), partelmat, nd*i, nd*j);
1546  }
1547  }
1548  else
1549  {
1550  if (Q)
1551  {
1552  norm *= Q->Eval(Trans, ip);
1553  }
1554  partelmat *= norm;
1555  for (int k = 0; k < vdim; k++)
1556  {
1557  elmat.AddMatrix(partelmat, nd*k, nd*k);
1558  }
1559  }
1560  }
1561 }
1562 
1563 void VectorMassIntegrator::AssembleElementMatrix2(
1564  const FiniteElement &trial_fe, const FiniteElement &test_fe,
1566 {
1567  int tr_nd = trial_fe.GetDof();
1568  int te_nd = test_fe.GetDof();
1569 
1570  double norm;
1571 
1572  // If vdim is not set, set it to the space dimension
1573  vdim = (vdim == -1) ? Trans.GetSpaceDim() : vdim;
1574 
1575  elmat.SetSize(te_nd*vdim, tr_nd*vdim);
1576  shape.SetSize(tr_nd);
1577  te_shape.SetSize(te_nd);
1578  partelmat.SetSize(te_nd, tr_nd);
1579  if (VQ)
1580  {
1581  vec.SetSize(vdim);
1582  }
1583  else if (MQ)
1584  {
1585  mcoeff.SetSize(vdim);
1586  }
1587 
1588  const IntegrationRule *ir = IntRule;
1589  if (ir == NULL)
1590  {
1591  int order = (trial_fe.GetOrder() + test_fe.GetOrder() +
1592  Trans.OrderW() + Q_order);
1593 
1594  if (trial_fe.Space() == FunctionSpace::rQk)
1595  {
1596  ir = &RefinedIntRules.Get(trial_fe.GetGeomType(), order);
1597  }
1598  else
1599  {
1600  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
1601  }
1602  }
1603 
1604  elmat = 0.0;
1605  for (int s = 0; s < ir->GetNPoints(); s++)
1606  {
1607  const IntegrationPoint &ip = ir->IntPoint(s);
1608  trial_fe.CalcShape(ip, shape);
1609  test_fe.CalcShape(ip, te_shape);
1610 
1611  Trans.SetIntPoint(&ip);
1612  norm = ip.weight * Trans.Weight();
1613 
1614  MultVWt(te_shape, shape, partelmat);
1615 
1616  if (VQ)
1617  {
1618  VQ->Eval(vec, Trans, ip);
1619  for (int k = 0; k < vdim; k++)
1620  {
1621  elmat.AddMatrix(norm*vec(k), partelmat, te_nd*k, tr_nd*k);
1622  }
1623  }
1624  else if (MQ)
1625  {
1626  MQ->Eval(mcoeff, Trans, ip);
1627  for (int i = 0; i < vdim; i++)
1628  for (int j = 0; j < vdim; j++)
1629  {
1630  elmat.AddMatrix(norm*mcoeff(i,j), partelmat, te_nd*i, tr_nd*j);
1631  }
1632  }
1633  else
1634  {
1635  if (Q)
1636  {
1637  norm *= Q->Eval(Trans, ip);
1638  }
1639  partelmat *= norm;
1640  for (int k = 0; k < vdim; k++)
1641  {
1642  elmat.AddMatrix(partelmat, te_nd*k, tr_nd*k);
1643  }
1644  }
1645  }
1646 }
1647 
1648 void VectorFEDivergenceIntegrator::AssembleElementMatrix2(
1649  const FiniteElement &trial_fe, const FiniteElement &test_fe,
1651 {
1652  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
1653 
1654 #ifdef MFEM_THREAD_SAFE
1655  Vector divshape(trial_nd), shape(test_nd);
1656 #else
1657  divshape.SetSize(trial_nd);
1658  shape.SetSize(test_nd);
1659 #endif
1660 
1661  elmat.SetSize(test_nd, trial_nd);
1662 
1663  const IntegrationRule *ir = IntRule;
1664  if (ir == NULL)
1665  {
1666  int order = trial_fe.GetOrder() + test_fe.GetOrder() - 1; // <--
1667  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
1668  }
1669 
1670  elmat = 0.0;
1671  for (i = 0; i < ir->GetNPoints(); i++)
1672  {
1673  const IntegrationPoint &ip = ir->IntPoint(i);
1674  trial_fe.CalcDivShape(ip, divshape);
1675  Trans.SetIntPoint(&ip);
1676  test_fe.CalcPhysShape(Trans, shape);
1677  double w = ip.weight;
1678  if (Q)
1679  {
1680  Trans.SetIntPoint(&ip);
1681  w *= Q->Eval(Trans, ip);
1682  }
1683  shape *= w;
1684  AddMultVWt(shape, divshape, elmat);
1685  }
1686 }
1687 
1688 void VectorFEWeakDivergenceIntegrator::AssembleElementMatrix2(
1689  const FiniteElement &trial_fe, const FiniteElement &test_fe,
1691 {
1692  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
1693  int dim = trial_fe.GetDim();
1694 
1695  MFEM_ASSERT(test_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
1696  test_fe.GetMapType() == mfem::FiniteElement::VALUE &&
1698  "Trial space must be H(Curl) and test space must be H_1");
1699 
1700 #ifdef MFEM_THREAD_SAFE
1701  DenseMatrix dshape(test_nd, dim);
1702  DenseMatrix dshapedxt(test_nd, dim);
1703  DenseMatrix vshape(trial_nd, dim);
1704  DenseMatrix invdfdx(dim);
1705 #else
1706  dshape.SetSize(test_nd, dim);
1707  dshapedxt.SetSize(test_nd, dim);
1708  vshape.SetSize(trial_nd, dim);
1709  invdfdx.SetSize(dim);
1710 #endif
1711 
1712  elmat.SetSize(test_nd, trial_nd);
1713 
1714  const IntegrationRule *ir = IntRule;
1715  if (ir == NULL)
1716  {
1717  // The integrand on the reference element is:
1718  // -( Q/det(J) ) u_hat^T adj(J) adj(J)^T grad_hat(v_hat).
1719  //
1720  // For Trans in (P_k)^d, v_hat in P_l, u_hat in ND_m, and dim=sdim=d>=1
1721  // - J_{ij} is in P_{k-1}, so adj(J)_{ij} is in P_{(d-1)*(k-1)}
1722  // - so adj(J)^T grad_hat(v_hat) is in (P_{(d-1)*(k-1)+(l-1)})^d
1723  // - u_hat is in (P_m)^d
1724  // - adj(J)^T u_hat is in (P_{(d-1)*(k-1)+m})^d
1725  // - and u_hat^T adj(J) adj(J)^T grad_hat(v_hat) is in P_n with
1726  // n = 2*(d-1)*(k-1)+(l-1)+m
1727  //
1728  // For Trans in (Q_k)^d, v_hat in Q_l, u_hat in ND_m, and dim=sdim=d>1
1729  // - J_{i*}, J's i-th row, is in ( Q_{k-1,k,k}, Q_{k,k-1,k}, Q_{k,k,k-1} )
1730  // - adj(J)_{*j} is in ( Q_{s,s-1,s-1}, Q_{s-1,s,s-1}, Q_{s-1,s-1,s} )
1731  // with s = (d-1)*k
1732  // - adj(J)^T grad_hat(v_hat) is in Q_{(d-1)*k+(l-1)}
1733  // - u_hat is in ( Q_{m-1,m,m}, Q_{m,m-1,m}, Q_{m,m,m-1} )
1734  // - adj(J)^T u_hat is in Q_{(d-1)*k+(m-1)}
1735  // - and u_hat^T adj(J) adj(J)^T grad_hat(v_hat) is in Q_n with
1736  // n = 2*(d-1)*k+(l-1)+(m-1)
1737  //
1738  // In the next formula we use the expressions for n with k=1, which means
1739  // that the term Q/det(J) is disregarded:
1740  int ir_order = (trial_fe.Space() == FunctionSpace::Pk) ?
1741  (trial_fe.GetOrder() + test_fe.GetOrder() - 1) :
1742  (trial_fe.GetOrder() + test_fe.GetOrder() + 2*(dim-2));
1743  ir = &IntRules.Get(trial_fe.GetGeomType(), ir_order);
1744  }
1745 
1746  elmat = 0.0;
1747  for (i = 0; i < ir->GetNPoints(); i++)
1748  {
1749  const IntegrationPoint &ip = ir->IntPoint(i);
1750  test_fe.CalcDShape(ip, dshape);
1751 
1752  Trans.SetIntPoint(&ip);
1753  CalcAdjugate(Trans.Jacobian(), invdfdx);
1754  Mult(dshape, invdfdx, dshapedxt);
1755 
1756  trial_fe.CalcVShape(Trans, vshape);
1757 
1758  double w = ip.weight;
1759 
1760  if (Q)
1761  {
1762  w *= Q->Eval(Trans, ip);
1763  }
1764  dshapedxt *= -w;
1765 
1766  AddMultABt(dshapedxt, vshape, elmat);
1767  }
1768 }
1769 
1770 void VectorFECurlIntegrator::AssembleElementMatrix2(
1771  const FiniteElement &trial_fe, const FiniteElement &test_fe,
1773 {
1774  int trial_nd = trial_fe.GetDof(), test_nd = test_fe.GetDof(), i;
1775  int dim = trial_fe.GetDim();
1776  int dimc = (dim == 3) ? 3 : 1;
1777 
1778  MFEM_ASSERT(trial_fe.GetMapType() == mfem::FiniteElement::H_CURL ||
1780  "At least one of the finite elements must be in H(Curl)");
1781 
1782  int curl_nd, vec_nd;
1783  if ( trial_fe.GetMapType() == mfem::FiniteElement::H_CURL )
1784  {
1785  curl_nd = trial_nd;
1786  vec_nd = test_nd;
1787  }
1788  else
1789  {
1790  curl_nd = test_nd;
1791  vec_nd = trial_nd;
1792  }
1793 
1794 #ifdef MFEM_THREAD_SAFE
1795  DenseMatrix curlshapeTrial(curl_nd, dimc);
1796  DenseMatrix curlshapeTrial_dFT(curl_nd, dimc);
1797  DenseMatrix vshapeTest(vec_nd, dimc);
1798 #else
1799  curlshapeTrial.SetSize(curl_nd, dimc);
1800  curlshapeTrial_dFT.SetSize(curl_nd, dimc);
1801  vshapeTest.SetSize(vec_nd, dimc);
1802 #endif
1803  Vector shapeTest(vshapeTest.GetData(), vec_nd);
1804 
1805  elmat.SetSize(test_nd, trial_nd);
1806 
1807  const IntegrationRule *ir = IntRule;
1808  if (ir == NULL)
1809  {
1810  int order = trial_fe.GetOrder() + test_fe.GetOrder() - 1; // <--
1811  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
1812  }
1813 
1814  elmat = 0.0;
1815  for (i = 0; i < ir->GetNPoints(); i++)
1816  {
1817  const IntegrationPoint &ip = ir->IntPoint(i);
1818 
1819  Trans.SetIntPoint(&ip);
1820  if (dim == 3)
1821  {
1822  if ( trial_fe.GetMapType() == mfem::FiniteElement::H_CURL )
1823  {
1824  trial_fe.CalcCurlShape(ip, curlshapeTrial);
1825  test_fe.CalcVShape(Trans, vshapeTest);
1826  }
1827  else
1828  {
1829  test_fe.CalcCurlShape(ip, curlshapeTrial);
1830  trial_fe.CalcVShape(Trans, vshapeTest);
1831  }
1832  MultABt(curlshapeTrial, Trans.Jacobian(), curlshapeTrial_dFT);
1833  }
1834  else
1835  {
1836  if ( trial_fe.GetMapType() == mfem::FiniteElement::H_CURL )
1837  {
1838  trial_fe.CalcCurlShape(ip, curlshapeTrial_dFT);
1839  test_fe.CalcShape(ip, shapeTest);
1840  }
1841  else
1842  {
1843  test_fe.CalcCurlShape(ip, curlshapeTrial_dFT);
1844  trial_fe.CalcShape(ip, shapeTest);
1845  }
1846  }
1847 
1848  double w = ip.weight;
1849 
1850  if (Q)
1851  {
1852  w *= Q->Eval(Trans, ip);
1853  }
1854  // Note: shapeTest points to the same data as vshapeTest
1855  vshapeTest *= w;
1856  if ( trial_fe.GetMapType() == mfem::FiniteElement::H_CURL )
1857  {
1858  AddMultABt(vshapeTest, curlshapeTrial_dFT, elmat);
1859  }
1860  else
1861  {
1862  AddMultABt(curlshapeTrial_dFT, vshapeTest, elmat);
1863  }
1864  }
1865 }
1866 
1867 void DerivativeIntegrator::AssembleElementMatrix2 (
1868  const FiniteElement &trial_fe,
1869  const FiniteElement &test_fe,
1871  DenseMatrix &elmat)
1872 {
1873  int dim = trial_fe.GetDim();
1874  int trial_nd = trial_fe.GetDof();
1875  int test_nd = test_fe.GetDof();
1876  int spaceDim = Trans.GetSpaceDim();
1877 
1878  int i, l;
1879  double det;
1880 
1881  elmat.SetSize (test_nd,trial_nd);
1882  dshape.SetSize (trial_nd,dim);
1883  dshapedxt.SetSize(trial_nd, spaceDim);
1884  dshapedxi.SetSize(trial_nd);
1885  invdfdx.SetSize(dim, spaceDim);
1886  shape.SetSize (test_nd);
1887 
1888  const IntegrationRule *ir = IntRule;
1889  if (ir == NULL)
1890  {
1891  int order;
1892  if (trial_fe.Space() == FunctionSpace::Pk)
1893  {
1894  order = trial_fe.GetOrder() + test_fe.GetOrder() - 1;
1895  }
1896  else
1897  {
1898  order = trial_fe.GetOrder() + test_fe.GetOrder() + dim;
1899  }
1900 
1901  if (trial_fe.Space() == FunctionSpace::rQk)
1902  {
1903  ir = &RefinedIntRules.Get(trial_fe.GetGeomType(), order);
1904  }
1905  else
1906  {
1907  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
1908  }
1909  }
1910 
1911  elmat = 0.0;
1912  for (i = 0; i < ir->GetNPoints(); i++)
1913  {
1914  const IntegrationPoint &ip = ir->IntPoint(i);
1915 
1916  trial_fe.CalcDShape(ip, dshape);
1917 
1918  Trans.SetIntPoint (&ip);
1919  CalcInverse (Trans.Jacobian(), invdfdx);
1920  det = Trans.Weight();
1921  Mult (dshape, invdfdx, dshapedxt);
1922 
1923  test_fe.CalcShape(ip, shape);
1924 
1925  for (l = 0; l < trial_nd; l++)
1926  {
1927  dshapedxi(l) = dshapedxt(l,xi);
1928  }
1929 
1930  shape *= Q->Eval(Trans,ip) * det * ip.weight;
1931  AddMultVWt (shape, dshapedxi, elmat);
1932  }
1933 }
1934 
1935 void CurlCurlIntegrator::AssembleElementMatrix
1937  DenseMatrix &elmat )
1938 {
1939  int nd = el.GetDof();
1940  dim = el.GetDim();
1941  int dimc = el.GetCurlDim();
1942  double w;
1943 
1944 #ifdef MFEM_THREAD_SAFE
1945  Vector D;
1946  DenseMatrix curlshape(nd,dimc), curlshape_dFt(nd,dimc), M;
1947 #else
1948  curlshape.SetSize(nd,dimc);
1949  curlshape_dFt.SetSize(nd,dimc);
1950 #endif
1951  elmat.SetSize(nd);
1952  if (MQ) { M.SetSize(dimc); }
1953  if (DQ) { D.SetSize(dimc); }
1954 
1955  const IntegrationRule *ir = IntRule;
1956  if (ir == NULL)
1957  {
1958  int order;
1959  if (el.Space() == FunctionSpace::Pk)
1960  {
1961  order = 2*el.GetOrder() - 2;
1962  }
1963  else
1964  {
1965  order = 2*el.GetOrder();
1966  }
1967 
1968  ir = &IntRules.Get(el.GetGeomType(), order);
1969  }
1970 
1971  elmat = 0.0;
1972  for (int i = 0; i < ir->GetNPoints(); i++)
1973  {
1974  const IntegrationPoint &ip = ir->IntPoint(i);
1975 
1976  Trans.SetIntPoint (&ip);
1977 
1978  w = ip.weight * Trans.Weight();
1979  el.CalcPhysCurlShape(Trans, curlshape_dFt);
1980 
1981  if (MQ)
1982  {
1983  MQ->Eval(M, Trans, ip);
1984  M *= w;
1985  Mult(curlshape_dFt, M, curlshape);
1986  AddMultABt(curlshape, curlshape_dFt, elmat);
1987  }
1988  else if (DQ)
1989  {
1990  DQ->Eval(D, Trans, ip);
1991  D *= w;
1992  AddMultADAt(curlshape_dFt, D, elmat);
1993  }
1994  else if (Q)
1995  {
1996  w *= Q->Eval(Trans, ip);
1997  AddMult_a_AAt(w, curlshape_dFt, elmat);
1998  }
1999  else
2000  {
2001  AddMult_a_AAt(w, curlshape_dFt, elmat);
2002  }
2003  }
2004 }
2005 
2006 void CurlCurlIntegrator::AssembleElementMatrix2(const FiniteElement &trial_fe,
2007  const FiniteElement &test_fe,
2009  DenseMatrix &elmat)
2010 {
2011  int tr_nd = trial_fe.GetDof();
2012  int te_nd = test_fe.GetDof();
2013  dim = trial_fe.GetDim();
2014  int dimc = trial_fe.GetCurlDim();
2015  double w;
2016 
2017 #ifdef MFEM_THREAD_SAFE
2018  Vector D;
2019  DenseMatrix curlshape(tr_nd,dimc), curlshape_dFt(tr_nd,dimc), M;
2020  DenseMatrix te_curlshape(te_nd,dimc), te_curlshape_dFt(te_nd,dimc);
2021 #else
2022  curlshape.SetSize(tr_nd,dimc);
2023  curlshape_dFt.SetSize(tr_nd,dimc);
2024  te_curlshape.SetSize(te_nd,dimc);
2025  te_curlshape_dFt.SetSize(te_nd,dimc);
2026 #endif
2027  elmat.SetSize(te_nd, tr_nd);
2028 
2029  if (MQ) { M.SetSize(dimc); }
2030  if (DQ) { D.SetSize(dimc); }
2031 
2032  const IntegrationRule *ir = IntRule;
2033  if (ir == NULL)
2034  {
2035  int order;
2036  if (trial_fe.Space() == FunctionSpace::Pk)
2037  {
2038  order = test_fe.GetOrder() + trial_fe.GetOrder() - 2;
2039  }
2040  else
2041  {
2042  order = test_fe.GetOrder() + trial_fe.GetOrder() + trial_fe.GetDim() - 1;
2043  }
2044  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
2045  }
2046 
2047  elmat = 0.0;
2048  for (int i = 0; i < ir->GetNPoints(); i++)
2049  {
2050  const IntegrationPoint &ip = ir->IntPoint(i);
2051 
2052  Trans.SetIntPoint(&ip);
2053 
2054  w = ip.weight * Trans.Weight();
2055  trial_fe.CalcPhysCurlShape(Trans, curlshape_dFt);
2056  test_fe.CalcPhysCurlShape(Trans, te_curlshape_dFt);
2057 
2058  if (MQ)
2059  {
2060  MQ->Eval(M, Trans, ip);
2061  M *= w;
2062  Mult(te_curlshape_dFt, M, te_curlshape);
2063  AddMultABt(te_curlshape, curlshape_dFt, elmat);
2064  }
2065  else if (DQ)
2066  {
2067  DQ->Eval(D, Trans, ip);
2068  D *= w;
2069  AddMultADBt(te_curlshape_dFt,D,curlshape_dFt,elmat);
2070  }
2071  else
2072  {
2073  if (Q)
2074  {
2075  w *= Q->Eval(Trans, ip);
2076  }
2077  curlshape_dFt *= w;
2078  AddMultABt(te_curlshape_dFt, curlshape_dFt, elmat);
2079  }
2080  }
2081 }
2082 
2083 void CurlCurlIntegrator
2084 ::ComputeElementFlux(const FiniteElement &el, ElementTransformation &Trans,
2085  Vector &u, const FiniteElement &fluxelem, Vector &flux,
2086  bool with_coef, const IntegrationRule *ir)
2087 {
2088 #ifdef MFEM_THREAD_SAFE
2089  DenseMatrix projcurl;
2090 #endif
2091 
2092  MFEM_VERIFY(ir == NULL, "Integration rule (ir) must be NULL")
2093 
2094  fluxelem.ProjectCurl(el, Trans, projcurl);
2095 
2096  flux.SetSize(projcurl.Height());
2097  projcurl.Mult(u, flux);
2098 
2099  // TODO: Q, wcoef?
2100 }
2101 
2102 double CurlCurlIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
2104  Vector &flux, Vector *d_energy)
2105 {
2106  int nd = fluxelem.GetDof();
2107  dim = fluxelem.GetDim();
2108 
2109 #ifdef MFEM_THREAD_SAFE
2110  DenseMatrix vshape;
2111 #endif
2112  vshape.SetSize(nd, dim);
2113  pointflux.SetSize(dim);
2114  if (d_energy) { vec.SetSize(dim); }
2115 
2116  int order = 2 * fluxelem.GetOrder(); // <--
2117  const IntegrationRule &ir = IntRules.Get(fluxelem.GetGeomType(), order);
2118 
2119  double energy = 0.0;
2120  if (d_energy) { *d_energy = 0.0; }
2121 
2122  Vector* pfluxes = NULL;
2123  if (d_energy)
2124  {
2125  pfluxes = new Vector[ir.GetNPoints()];
2126  }
2127 
2128  for (int i = 0; i < ir.GetNPoints(); i++)
2129  {
2130  const IntegrationPoint &ip = ir.IntPoint(i);
2131  Trans.SetIntPoint(&ip);
2132 
2133  fluxelem.CalcVShape(Trans, vshape);
2134  // fluxelem.CalcVShape(ip, vshape);
2135  vshape.MultTranspose(flux, pointflux);
2136 
2137  double w = Trans.Weight() * ip.weight;
2138 
2139  double e = w * (pointflux * pointflux);
2140 
2141  if (Q)
2142  {
2143  // TODO
2144  }
2145 
2146  energy += e;
2147 
2148 #if ANISO_EXPERIMENTAL
2149  if (d_energy)
2150  {
2151  pfluxes[i].SetSize(dim);
2152  Trans.Jacobian().MultTranspose(pointflux, pfluxes[i]);
2153 
2154  /*
2155  DenseMatrix Jadj(dim, dim);
2156  CalcAdjugate(Trans.Jacobian(), Jadj);
2157  pfluxes[i].SetSize(dim);
2158  Jadj.Mult(pointflux, pfluxes[i]);
2159  */
2160 
2161  // pfluxes[i] = pointflux;
2162  }
2163 #endif
2164  }
2165 
2166  if (d_energy)
2167  {
2168 #if ANISO_EXPERIMENTAL
2169  *d_energy = 0.0;
2170  Vector tmp;
2171 
2172  int n = (int) round(pow(ir.GetNPoints(), 1.0/3.0));
2173  MFEM_ASSERT(n*n*n == ir.GetNPoints(), "");
2174 
2175  // hack: get total variation of 'pointflux' in the x,y,z directions
2176  for (int k = 0; k < n; k++)
2177  for (int l = 0; l < n; l++)
2178  for (int m = 0; m < n; m++)
2179  {
2180  Vector &vec = pfluxes[(k*n + l)*n + m];
2181  if (m > 0)
2182  {
2183  tmp = vec; tmp -= pfluxes[(k*n + l)*n + (m-1)];
2184  (*d_energy)[0] += (tmp * tmp);
2185  }
2186  if (l > 0)
2187  {
2188  tmp = vec; tmp -= pfluxes[(k*n + (l-1))*n + m];
2189  (*d_energy)[1] += (tmp * tmp);
2190  }
2191  if (k > 0)
2192  {
2193  tmp = vec; tmp -= pfluxes[((k-1)*n + l)*n + m];
2194  (*d_energy)[2] += (tmp * tmp);
2195  }
2196  }
2197 #else
2198  *d_energy = 1.0;
2199 #endif
2200 
2201  delete [] pfluxes;
2202  }
2203 
2204  return energy;
2205 }
2206 
2207 void VectorCurlCurlIntegrator::AssembleElementMatrix(
2209 {
2210  int dim = el.GetDim();
2211  int dof = el.GetDof();
2212  int cld = (dim*(dim-1))/2;
2213 
2214 #ifdef MFEM_THREAD_SAFE
2215  DenseMatrix dshape_hat(dof, dim), dshape(dof, dim);
2216  DenseMatrix curlshape(dim*dof, cld), Jadj(dim);
2217 #else
2218  dshape_hat.SetSize(dof, dim);
2219  dshape.SetSize(dof, dim);
2220  curlshape.SetSize(dim*dof, cld);
2221  Jadj.SetSize(dim);
2222 #endif
2223 
2224  const IntegrationRule *ir = IntRule;
2225  if (ir == NULL)
2226  {
2227  // use the same integration rule as diffusion
2228  int order = 2 * Trans.OrderGrad(&el);
2229  ir = &IntRules.Get(el.GetGeomType(), order);
2230  }
2231 
2232  elmat.SetSize(dof*dim);
2233  elmat = 0.0;
2234  for (int i = 0; i < ir->GetNPoints(); i++)
2235  {
2236  const IntegrationPoint &ip = ir->IntPoint(i);
2237  el.CalcDShape(ip, dshape_hat);
2238 
2239  Trans.SetIntPoint(&ip);
2240  CalcAdjugate(Trans.Jacobian(), Jadj);
2241  double w = ip.weight / Trans.Weight();
2242 
2243  Mult(dshape_hat, Jadj, dshape);
2244  dshape.GradToCurl(curlshape);
2245 
2246  if (Q)
2247  {
2248  w *= Q->Eval(Trans, ip);
2249  }
2250 
2251  AddMult_a_AAt(w, curlshape, elmat);
2252  }
2253 }
2254 
2255 double VectorCurlCurlIntegrator::GetElementEnergy(
2256  const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun)
2257 {
2258  int dim = el.GetDim();
2259  int dof = el.GetDof();
2260 
2261 #ifdef MFEM_THREAD_SAFE
2262  DenseMatrix dshape_hat(dof, dim), Jadj(dim), grad_hat(dim), grad(dim);
2263 #else
2264  dshape_hat.SetSize(dof, dim);
2265 
2266  Jadj.SetSize(dim);
2267  grad_hat.SetSize(dim);
2268  grad.SetSize(dim);
2269 #endif
2270  DenseMatrix elfun_mat(elfun.GetData(), dof, dim);
2271 
2272  const IntegrationRule *ir = IntRule;
2273  if (ir == NULL)
2274  {
2275  // use the same integration rule as diffusion
2276  int order = 2 * Tr.OrderGrad(&el);
2277  ir = &IntRules.Get(el.GetGeomType(), order);
2278  }
2279 
2280  double energy = 0.;
2281  for (int i = 0; i < ir->GetNPoints(); i++)
2282  {
2283  const IntegrationPoint &ip = ir->IntPoint(i);
2284  el.CalcDShape(ip, dshape_hat);
2285 
2286  MultAtB(elfun_mat, dshape_hat, grad_hat);
2287 
2288  Tr.SetIntPoint(&ip);
2289  CalcAdjugate(Tr.Jacobian(), Jadj);
2290  double w = ip.weight / Tr.Weight();
2291 
2292  Mult(grad_hat, Jadj, grad);
2293 
2294  if (dim == 2)
2295  {
2296  double curl = grad(0,1) - grad(1,0);
2297  w *= curl * curl;
2298  }
2299  else
2300  {
2301  double curl_x = grad(2,1) - grad(1,2);
2302  double curl_y = grad(0,2) - grad(2,0);
2303  double curl_z = grad(1,0) - grad(0,1);
2304  w *= curl_x * curl_x + curl_y * curl_y + curl_z * curl_z;
2305  }
2306 
2307  if (Q)
2308  {
2309  w *= Q->Eval(Tr, ip);
2310  }
2311 
2312  energy += w;
2313  }
2314 
2315  elfun_mat.ClearExternalData();
2316 
2317  return 0.5 * energy;
2318 }
2319 
2320 void MixedCurlIntegrator::AssembleElementMatrix2(
2321  const FiniteElement &trial_fe, const FiniteElement &test_fe,
2323 {
2324  int dim = trial_fe.GetDim();
2325  int trial_dof = trial_fe.GetDof();
2326  int test_dof = test_fe.GetDof();
2327  int dimc = (dim == 3) ? 3 : 1;
2328 
2329  MFEM_VERIFY(trial_fe.GetMapType() == mfem::FiniteElement::H_CURL ||
2330  (dim == 2 && trial_fe.GetMapType() == mfem::FiniteElement::VALUE),
2331  "Trial finite element must be either 2D/3D H(Curl) or 2D H1");
2332  MFEM_VERIFY(test_fe.GetMapType() == mfem::FiniteElement::VALUE ||
2334  "Test finite element must be in H1/L2");
2335 
2336  bool spaceH1 = (trial_fe.GetMapType() == mfem::FiniteElement::VALUE);
2337 
2338  if (spaceH1)
2339  {
2340  dshape.SetSize(trial_dof,dim);
2341  curlshape.SetSize(trial_dof,dim);
2342  dimc = dim;
2343  }
2344  else
2345  {
2346  curlshape.SetSize(trial_dof,dimc);
2347  elmat_comp.SetSize(test_dof, trial_dof);
2348  }
2349  elmat.SetSize(dimc * test_dof, trial_dof);
2350  shape.SetSize(test_dof);
2351  elmat = 0.0;
2352 
2353  double c;
2354  Vector d_col;
2355  const IntegrationRule *ir = IntRule;
2356 
2357  if (ir == NULL)
2358  {
2359  int order = trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderJ();
2360  ir = &IntRules.Get(trial_fe.GetGeomType(), order);
2361  }
2362 
2363  for (int i = 0; i < ir->GetNPoints(); i++)
2364  {
2365  const IntegrationPoint &ip = ir->IntPoint(i);
2366  Trans.SetIntPoint(&ip);
2367  if (spaceH1)
2368  {
2369  trial_fe.CalcPhysDShape(Trans, dshape);
2370  dshape.GradToVectorCurl2D(curlshape);
2371  }
2372  else
2373  {
2374  trial_fe.CalcPhysCurlShape(Trans, curlshape);
2375  }
2376  test_fe.CalcPhysShape(Trans, shape);
2377  c = ip.weight*Trans.Weight();
2378  if (Q)
2379  {
2380  c *= Q->Eval(Trans, ip);
2381  }
2382  shape *= c;
2383 
2384  for (int d = 0; d < dimc; ++d)
2385  {
2386  double * curldata = &(curlshape.GetData())[d*trial_dof];
2387  for (int jj = 0; jj < trial_dof; ++jj)
2388  {
2389  for (int ii = 0; ii < test_dof; ++ii)
2390  {
2391  elmat(d * test_dof + ii, jj) += shape(ii) * curldata[jj];
2392  }
2393  }
2394  }
2395  }
2396 }
2397 
2398 
2399 void VectorFEMassIntegrator::AssembleElementMatrix(
2400  const FiniteElement &el,
2402  DenseMatrix &elmat)
2403 {
2404  int dof = el.GetDof();
2405  int spaceDim = Trans.GetSpaceDim();
2406  int vdim = std::max(spaceDim, el.GetVDim());
2407 
2408  double w;
2409 
2410 #ifdef MFEM_THREAD_SAFE
2411  Vector D(DQ ? DQ->GetVDim() : 0);
2412  DenseMatrix trial_vshape(dof, vdim);
2413  DenseMatrix K(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
2414 #else
2415  trial_vshape.SetSize(dof, vdim);
2416  D.SetSize(DQ ? DQ->GetVDim() : 0);
2417  K.SetSize(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
2418 #endif
2419  DenseMatrix tmp(trial_vshape.Height(), K.Width());
2420 
2421  elmat.SetSize(dof);
2422  elmat = 0.0;
2423 
2424  const IntegrationRule *ir = IntRule;
2425  if (ir == NULL)
2426  {
2427  // int order = 2 * el.GetOrder();
2428  int order = Trans.OrderW() + 2 * el.GetOrder();
2429  ir = &IntRules.Get(el.GetGeomType(), order);
2430  }
2431 
2432  for (int i = 0; i < ir->GetNPoints(); i++)
2433  {
2434  const IntegrationPoint &ip = ir->IntPoint(i);
2435 
2436  Trans.SetIntPoint (&ip);
2437 
2438  el.CalcVShape(Trans, trial_vshape);
2439 
2440  w = ip.weight * Trans.Weight();
2441  if (MQ)
2442  {
2443  MQ->Eval(K, Trans, ip);
2444  K *= w;
2445  Mult(trial_vshape,K,tmp);
2446  AddMultABt(tmp,trial_vshape,elmat);
2447  }
2448  else if (DQ)
2449  {
2450  DQ->Eval(D, Trans, ip);
2451  D *= w;
2452  AddMultADAt(trial_vshape, D, elmat);
2453  }
2454  else
2455  {
2456  if (Q)
2457  {
2458  w *= Q -> Eval (Trans, ip);
2459  }
2460  AddMult_a_AAt (w, trial_vshape, elmat);
2461  }
2462  }
2463 }
2464 
2465 void VectorFEMassIntegrator::AssembleElementMatrix2(
2466  const FiniteElement &trial_fe, const FiniteElement &test_fe,
2468 {
2469  if (test_fe.GetRangeType() == FiniteElement::SCALAR
2470  && trial_fe.GetRangeType() == FiniteElement::VECTOR)
2471  {
2472  // assume test_fe is scalar FE and trial_fe is vector FE
2473  int spaceDim = Trans.GetSpaceDim();
2474  int vdim = std::max(spaceDim, trial_fe.GetVDim());
2475  int trial_dof = trial_fe.GetDof();
2476  int test_dof = test_fe.GetDof();
2477  double w;
2478 
2479 #ifdef MFEM_THREAD_SAFE
2480  DenseMatrix trial_vshape(trial_dof, spaceDim);
2481  Vector shape(test_dof);
2482  Vector D(DQ ? DQ->GetVDim() : 0);
2483  DenseMatrix K(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
2484 #else
2485  trial_vshape.SetSize(trial_dof, spaceDim);
2486  shape.SetSize(test_dof);
2487  D.SetSize(DQ ? DQ->GetVDim() : 0);
2488  K.SetSize(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
2489 #endif
2490 
2491  elmat.SetSize(vdim*test_dof, trial_dof);
2492 
2493  const IntegrationRule *ir = IntRule;
2494  if (ir == NULL)
2495  {
2496  int order = (Trans.OrderW() + test_fe.GetOrder() + trial_fe.GetOrder());
2497  ir = &IntRules.Get(test_fe.GetGeomType(), order);
2498  }
2499 
2500  elmat = 0.0;
2501  for (int i = 0; i < ir->GetNPoints(); i++)
2502  {
2503  const IntegrationPoint &ip = ir->IntPoint(i);
2504 
2505  Trans.SetIntPoint (&ip);
2506 
2507  trial_fe.CalcVShape(Trans, trial_vshape);
2508  test_fe.CalcShape(ip, shape);
2509 
2510  w = ip.weight * Trans.Weight();
2511  if (DQ)
2512  {
2513  DQ->Eval(D, Trans, ip);
2514  D *= w;
2515  for (int d = 0; d < vdim; d++)
2516  {
2517  for (int j = 0; j < test_dof; j++)
2518  {
2519  for (int k = 0; k < trial_dof; k++)
2520  {
2521  elmat(d * test_dof + j, k) +=
2522  shape(j) * D(d) * trial_vshape(k, d);
2523  }
2524  }
2525  }
2526  }
2527  else if (MQ)
2528  {
2529  MQ->Eval(K, Trans, ip);
2530  K *= w;
2531  for (int d = 0; d < vdim; d++)
2532  {
2533  for (int j = 0; j < test_dof; j++)
2534  {
2535  for (int k = 0; k < trial_dof; k++)
2536  {
2537  double Kv = 0.0;
2538  for (int vd = 0; vd < spaceDim; vd++)
2539  {
2540  Kv += K(d, vd) * trial_vshape(k, vd);
2541  }
2542  elmat(d * test_dof + j, k) += shape(j) * Kv;
2543  }
2544  }
2545  }
2546  }
2547  else
2548  {
2549  if (Q)
2550  {
2551  w *= Q->Eval(Trans, ip);
2552  }
2553  for (int d = 0; d < vdim; d++)
2554  {
2555  for (int j = 0; j < test_dof; j++)
2556  {
2557  for (int k = 0; k < trial_dof; k++)
2558  {
2559  elmat(d * test_dof + j, k) +=
2560  w * shape(j) * trial_vshape(k, d);
2561  }
2562  }
2563  }
2564  }
2565  }
2566  }
2567  else if (test_fe.GetRangeType() == FiniteElement::VECTOR
2568  && trial_fe.GetRangeType() == FiniteElement::VECTOR)
2569  {
2570  // assume both test_fe and trial_fe are vector FE
2571  int spaceDim = Trans.GetSpaceDim();
2572  int trial_vdim = std::max(spaceDim, trial_fe.GetVDim());
2573  int test_vdim = std::max(spaceDim, test_fe.GetVDim());
2574  int trial_dof = trial_fe.GetDof();
2575  int test_dof = test_fe.GetDof();
2576  double w;
2577 
2578 #ifdef MFEM_THREAD_SAFE
2579  DenseMatrix trial_vshape(trial_dof,trial_vdim);
2580  DenseMatrix test_vshape(test_dof,test_vdim);
2581  Vector D(DQ ? DQ->GetVDim() : 0);
2582  DenseMatrix K(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
2583 #else
2584  trial_vshape.SetSize(trial_dof,trial_vdim);
2585  test_vshape.SetSize(test_dof,test_vdim);
2586  D.SetSize(DQ ? DQ->GetVDim() : 0);
2587  K.SetSize(MQ ? MQ->GetVDim() : 0, MQ ? MQ->GetVDim() : 0);
2588 #endif
2589  DenseMatrix tmp(test_vshape.Height(), K.Width());
2590 
2591  elmat.SetSize (test_dof, trial_dof);
2592 
2593  const IntegrationRule *ir = IntRule;
2594  if (ir == NULL)
2595  {
2596  int order = (Trans.OrderW() + test_fe.GetOrder() + trial_fe.GetOrder());
2597  ir = &IntRules.Get(test_fe.GetGeomType(), order);
2598  }
2599 
2600  elmat = 0.0;
2601  for (int i = 0; i < ir->GetNPoints(); i++)
2602  {
2603  const IntegrationPoint &ip = ir->IntPoint(i);
2604 
2605  Trans.SetIntPoint (&ip);
2606 
2607  trial_fe.CalcVShape(Trans, trial_vshape);
2608  test_fe.CalcVShape(Trans, test_vshape);
2609 
2610  w = ip.weight * Trans.Weight();
2611  if (MQ)
2612  {
2613  MQ->Eval(K, Trans, ip);
2614  K *= w;
2615  Mult(test_vshape,K,tmp);
2616  AddMultABt(tmp,trial_vshape,elmat);
2617  }
2618  else if (DQ)
2619  {
2620  DQ->Eval(D, Trans, ip);
2621  D *= w;
2622  AddMultADBt(test_vshape,D,trial_vshape,elmat);
2623  }
2624  else
2625  {
2626  if (Q)
2627  {
2628  w *= Q -> Eval (Trans, ip);
2629  }
2630  AddMult_a_ABt(w,test_vshape,trial_vshape,elmat);
2631  }
2632  }
2633  }
2634  else
2635  {
2636  mfem_error("VectorFEMassIntegrator::AssembleElementMatrix2(...)\n"
2637  " is not implemented for given trial and test bases.");
2638  }
2639 }
2640 
2641 void VectorDivergenceIntegrator::AssembleElementMatrix2(
2642  const FiniteElement &trial_fe,
2643  const FiniteElement &test_fe,
2645  DenseMatrix &elmat)
2646 {
2647  dim = trial_fe.GetDim();
2648  int trial_dof = trial_fe.GetDof();
2649  int test_dof = test_fe.GetDof();
2650  double c;
2651 
2652  dshape.SetSize (trial_dof, dim);
2653  gshape.SetSize (trial_dof, dim);
2654  Jadj.SetSize (dim);
2655  divshape.SetSize (dim*trial_dof);
2656  shape.SetSize (test_dof);
2657 
2658  elmat.SetSize (test_dof, dim*trial_dof);
2659 
2660  const IntegrationRule *ir = IntRule ? IntRule : &GetRule(trial_fe, test_fe,
2661  Trans);
2662 
2663  elmat = 0.0;
2664 
2665  for (int i = 0; i < ir -> GetNPoints(); i++)
2666  {
2667  const IntegrationPoint &ip = ir->IntPoint(i);
2668 
2669  trial_fe.CalcDShape (ip, dshape);
2670  test_fe.CalcShape (ip, shape);
2671 
2672  Trans.SetIntPoint (&ip);
2673  CalcAdjugate(Trans.Jacobian(), Jadj);
2674 
2675  Mult (dshape, Jadj, gshape);
2676 
2677  gshape.GradToDiv (divshape);
2678 
2679  c = ip.weight;
2680  if (Q)
2681  {
2682  c *= Q -> Eval (Trans, ip);
2683  }
2684 
2685  // elmat += c * shape * divshape ^ t
2686  shape *= c;
2687  AddMultVWt (shape, divshape, elmat);
2688  }
2689 }
2690 
2692  const FiniteElement &trial_fe,
2693  const FiniteElement &test_fe,
2695 {
2696  int order = Trans.OrderGrad(&trial_fe) + test_fe.GetOrder() + Trans.OrderJ();
2697  return IntRules.Get(trial_fe.GetGeomType(), order);
2698 }
2699 
2700 
2701 void DivDivIntegrator::AssembleElementMatrix(
2702  const FiniteElement &el,
2704  DenseMatrix &elmat)
2705 {
2706  int dof = el.GetDof();
2707  double c;
2708 
2709 #ifdef MFEM_THREAD_SAFE
2710  Vector divshape(dof);
2711 #else
2712  divshape.SetSize(dof);
2713 #endif
2714  elmat.SetSize(dof);
2715 
2716  const IntegrationRule *ir = IntRule;
2717  if (ir == NULL)
2718  {
2719  int order = 2 * el.GetOrder() - 2; // <--- OK for RTk
2720  ir = &IntRules.Get(el.GetGeomType(), order);
2721  }
2722 
2723  elmat = 0.0;
2724 
2725  for (int i = 0; i < ir -> GetNPoints(); i++)
2726  {
2727  const IntegrationPoint &ip = ir->IntPoint(i);
2728 
2729  el.CalcDivShape (ip, divshape);
2730 
2731  Trans.SetIntPoint (&ip);
2732  c = ip.weight / Trans.Weight();
2733 
2734  if (Q)
2735  {
2736  c *= Q -> Eval (Trans, ip);
2737  }
2738 
2739  // elmat += c * divshape * divshape ^ t
2740  AddMult_a_VVt (c, divshape, elmat);
2741  }
2742 }
2743 
2744 void DivDivIntegrator::AssembleElementMatrix2(
2745  const FiniteElement &trial_fe,
2746  const FiniteElement &test_fe,
2748  DenseMatrix &elmat)
2749 {
2750  int tr_nd = trial_fe.GetDof();
2751  int te_nd = test_fe.GetDof();
2752  double c;
2753 
2754 #ifdef MFEM_THREAD_SAFE
2755  Vector divshape(tr_nd);
2756  Vector te_divshape(te_nd);
2757 #else
2758  divshape.SetSize(tr_nd);
2759  te_divshape.SetSize(te_nd);
2760 #endif
2761  elmat.SetSize(te_nd,tr_nd);
2762 
2763  const IntegrationRule *ir = IntRule;
2764  if (ir == NULL)
2765  {
2766  int order = 2 * max(test_fe.GetOrder(),
2767  trial_fe.GetOrder()) - 2; // <--- OK for RTk
2768  ir = &IntRules.Get(test_fe.GetGeomType(), order);
2769  }
2770 
2771  elmat = 0.0;
2772 
2773  for (int i = 0; i < ir -> GetNPoints(); i++)
2774  {
2775  const IntegrationPoint &ip = ir->IntPoint(i);
2776 
2777  trial_fe.CalcDivShape(ip,divshape);
2778  test_fe.CalcDivShape(ip,te_divshape);
2779 
2780  Trans.SetIntPoint (&ip);
2781  c = ip.weight / Trans.Weight();
2782 
2783  if (Q)
2784  {
2785  c *= Q -> Eval (Trans, ip);
2786  }
2787 
2788  te_divshape *= c;
2789  AddMultVWt(te_divshape, divshape, elmat);
2790  }
2791 }
2792 
2793 void VectorDiffusionIntegrator::AssembleElementMatrix(
2794  const FiniteElement &el,
2796  DenseMatrix &elmat)
2797 {
2798  const int dof = el.GetDof();
2799  dim = el.GetDim();
2800  sdim = Trans.GetSpaceDim();
2801 
2802  // If vdim is not set, set it to the space dimension;
2803  vdim = (vdim <= 0) ? sdim : vdim;
2804  const bool square = (dim == sdim);
2805 
2806  if (VQ)
2807  {
2808  vcoeff.SetSize(vdim);
2809  }
2810  else if (MQ)
2811  {
2812  mcoeff.SetSize(vdim);
2813  }
2814 
2815  dshape.SetSize(dof, dim);
2816  dshapedxt.SetSize(dof, sdim);
2817 
2818  elmat.SetSize(vdim * dof);
2819  pelmat.SetSize(dof);
2820 
2821  const IntegrationRule *ir = IntRule;
2822  if (ir == NULL)
2823  {
2824  ir = &DiffusionIntegrator::GetRule(el,el);
2825  }
2826 
2827  elmat = 0.0;
2828 
2829  for (int i = 0; i < ir -> GetNPoints(); i++)
2830  {
2831 
2832  const IntegrationPoint &ip = ir->IntPoint(i);
2833  el.CalcDShape(ip, dshape);
2834 
2835  Trans.SetIntPoint(&ip);
2836  double w = Trans.Weight();
2837  w = ip.weight / (square ? w : w*w*w);
2838  // AdjugateJacobian = / adj(J), if J is square
2839  // \ adj(J^t.J).J^t, otherwise
2840  Mult(dshape, Trans.AdjugateJacobian(), dshapedxt);
2841 
2842  if (VQ)
2843  {
2844  VQ->Eval(vcoeff, Trans, ip);
2845  for (int k = 0; k < vdim; ++k)
2846  {
2847  Mult_a_AAt(w*vcoeff(k), dshapedxt, pelmat);
2848  elmat.AddMatrix(pelmat, dof*k, dof*k);
2849  }
2850  }
2851  else if (MQ)
2852  {
2853  MQ->Eval(mcoeff, Trans, ip);
2854  for (int ii = 0; ii < vdim; ++ii)
2855  {
2856  for (int jj = 0; jj < vdim; ++jj)
2857  {
2858  Mult_a_AAt(w*mcoeff(ii,jj), dshapedxt, pelmat);
2859  elmat.AddMatrix(pelmat, dof*ii, dof*jj);
2860  }
2861  }
2862  }
2863  else
2864  {
2865  if (Q) { w *= Q->Eval(Trans, ip); }
2866  Mult_a_AAt(w, dshapedxt, pelmat);
2867  for (int k = 0; k < vdim; ++k)
2868  {
2869  elmat.AddMatrix(pelmat, dof*k, dof*k);
2870  }
2871  }
2872  }
2873 }
2874 
2875 void VectorDiffusionIntegrator::AssembleElementVector(
2876  const FiniteElement &el, ElementTransformation &Tr,
2877  const Vector &elfun, Vector &elvect)
2878 {
2879  const int dof = el.GetDof();
2880  dim = el.GetDim();
2881  sdim = Tr.GetSpaceDim();
2882 
2883  // If vdim is not set, set it to the space dimension;
2884  vdim = (vdim <= 0) ? sdim : vdim;
2885  const bool square = (dim == sdim);
2886 
2887  if (VQ)
2888  {
2889  vcoeff.SetSize(vdim);
2890  }
2891  else if (MQ)
2892  {
2893  mcoeff.SetSize(vdim);
2894  }
2895 
2896  dshape.SetSize(dof, dim);
2897  dshapedxt.SetSize(dof, dim);
2898  // pelmat.SetSize(dim);
2899 
2900  elvect.SetSize(dim*dof);
2901 
2902  // NOTE: DenseMatrix is in column-major order. This is consistent with
2903  // vectors ordered byNODES. In the resulting DenseMatrix, each column
2904  // corresponds to a particular vdim.
2905  DenseMatrix mat_in(elfun.GetData(), dof, dim);
2906  DenseMatrix mat_out(elvect.GetData(), dof, dim);
2907 
2908  const IntegrationRule *ir = IntRule;
2909  if (ir == NULL)
2910  {
2911  ir = &DiffusionIntegrator::GetRule(el,el);
2912  }
2913 
2914  elvect = 0.0;
2915  for (int i = 0; i < ir->GetNPoints(); i++)
2916  {
2917  const IntegrationPoint &ip = ir->IntPoint(i);
2918  el.CalcDShape(ip, dshape);
2919 
2920  Tr.SetIntPoint(&ip);
2921  double w = Tr.Weight();
2922  w = ip.weight / (square ? w : w*w*w);
2923  Mult(dshape, Tr.AdjugateJacobian(), dshapedxt);
2924  MultAAt(dshapedxt, pelmat);
2925 
2926  if (VQ)
2927  {
2928  VQ->Eval(vcoeff, Tr, ip);
2929  for (int k = 0; k < vdim; ++k)
2930  {
2931  pelmat *= w*vcoeff(k);
2932  const Vector vec_in(mat_in.GetColumn(k), dof);
2933  Vector vec_out(mat_out.GetColumn(k), dof);
2934  pelmat.AddMult(vec_in, vec_out);
2935  }
2936  }
2937  else if (MQ)
2938  {
2939  MQ->Eval(mcoeff, Tr, ip);
2940  for (int ii = 0; ii < vdim; ++ii)
2941  {
2942  Vector vec_out(mat_out.GetColumn(ii), dof);
2943  for (int jj = 0; jj < vdim; ++jj)
2944  {
2945  pelmat *= w*mcoeff(ii,jj);
2946  const Vector vec_in(mat_in.GetColumn(jj), dof);
2947  pelmat.Mult(vec_in, vec_out);
2948  }
2949  }
2950  }
2951  else
2952  {
2953  if (Q) { w *= Q->Eval(Tr, ip); }
2954  pelmat *= w;
2955  for (int k = 0; k < vdim; ++k)
2956  {
2957  const Vector vec_in(mat_in.GetColumn(k), dof);
2958  Vector vec_out(mat_out.GetColumn(k), dof);
2959  pelmat.AddMult(vec_in, vec_out);
2960  }
2961  }
2962  }
2963 }
2964 
2965 
2966 void ElasticityIntegrator::AssembleElementMatrix(
2968 {
2969  int dof = el.GetDof();
2970  int dim = el.GetDim();
2971  double w, L, M;
2972 
2973  MFEM_ASSERT(dim == Trans.GetSpaceDim(), "");
2974 
2975 #ifdef MFEM_THREAD_SAFE
2976  DenseMatrix dshape(dof, dim), gshape(dof, dim), pelmat(dof);
2977  Vector divshape(dim*dof);
2978 #else
2979  dshape.SetSize(dof, dim);
2980  gshape.SetSize(dof, dim);
2981  pelmat.SetSize(dof);
2982  divshape.SetSize(dim*dof);
2983 #endif
2984 
2985  elmat.SetSize(dof * dim);
2986 
2987  const IntegrationRule *ir = IntRule;
2988  if (ir == NULL)
2989  {
2990  int order = 2 * Trans.OrderGrad(&el); // correct order?
2991  ir = &IntRules.Get(el.GetGeomType(), order);
2992  }
2993 
2994  elmat = 0.0;
2995 
2996  for (int i = 0; i < ir -> GetNPoints(); i++)
2997  {
2998  const IntegrationPoint &ip = ir->IntPoint(i);
2999 
3000  el.CalcDShape(ip, dshape);
3001 
3002  Trans.SetIntPoint(&ip);
3003  w = ip.weight * Trans.Weight();
3004  Mult(dshape, Trans.InverseJacobian(), gshape);
3005  MultAAt(gshape, pelmat);
3006  gshape.GradToDiv (divshape);
3007 
3008  M = mu->Eval(Trans, ip);
3009  if (lambda)
3010  {
3011  L = lambda->Eval(Trans, ip);
3012  }
3013  else
3014  {
3015  L = q_lambda * M;
3016  M = q_mu * M;
3017  }
3018 
3019  if (L != 0.0)
3020  {
3021  AddMult_a_VVt(L * w, divshape, elmat);
3022  }
3023 
3024  if (M != 0.0)
3025  {
3026  for (int d = 0; d < dim; d++)
3027  {
3028  for (int k = 0; k < dof; k++)
3029  for (int l = 0; l < dof; l++)
3030  {
3031  elmat (dof*d+k, dof*d+l) += (M * w) * pelmat(k, l);
3032  }
3033  }
3034  for (int ii = 0; ii < dim; ii++)
3035  for (int jj = 0; jj < dim; jj++)
3036  {
3037  for (int kk = 0; kk < dof; kk++)
3038  for (int ll = 0; ll < dof; ll++)
3039  {
3040  elmat(dof*ii+kk, dof*jj+ll) +=
3041  (M * w) * gshape(kk, jj) * gshape(ll, ii);
3042  }
3043  }
3044  }
3045  }
3046 }
3047 
3048 void ElasticityIntegrator::ComputeElementFlux(
3050  Vector &u, const mfem::FiniteElement &fluxelem, Vector &flux,
3051  bool with_coef, const IntegrationRule *ir)
3052 {
3053  const int dof = el.GetDof();
3054  const int dim = el.GetDim();
3055  const int tdim = dim*(dim+1)/2; // num. entries in a symmetric tensor
3056  double L, M;
3057 
3058  MFEM_ASSERT(dim == 2 || dim == 3,
3059  "dimension is not supported: dim = " << dim);
3060  MFEM_ASSERT(dim == Trans.GetSpaceDim(), "");
3061  MFEM_ASSERT(fluxelem.GetMapType() == FiniteElement::VALUE, "");
3062  MFEM_ASSERT(dynamic_cast<const NodalFiniteElement*>(&fluxelem), "");
3063 
3064 #ifdef MFEM_THREAD_SAFE
3065  DenseMatrix dshape(dof, dim);
3066 #else
3067  dshape.SetSize(dof, dim);
3068 #endif
3069 
3070  double gh_data[9], grad_data[9];
3071  DenseMatrix gh(gh_data, dim, dim);
3072  DenseMatrix grad(grad_data, dim, dim);
3073 
3074  if (!ir)
3075  {
3076  ir = &fluxelem.GetNodes();
3077  }
3078  const int fnd = ir->GetNPoints();
3079  flux.SetSize(fnd * tdim);
3080 
3081  DenseMatrix loc_data_mat(u.GetData(), dof, dim);
3082  for (int i = 0; i < fnd; i++)
3083  {
3084  const IntegrationPoint &ip = ir->IntPoint(i);
3085  el.CalcDShape(ip, dshape);
3086  MultAtB(loc_data_mat, dshape, gh);
3087 
3088  Trans.SetIntPoint(&ip);
3089  Mult(gh, Trans.InverseJacobian(), grad);
3090 
3091  M = mu->Eval(Trans, ip);
3092  if (lambda)
3093  {
3094  L = lambda->Eval(Trans, ip);
3095  }
3096  else
3097  {
3098  L = q_lambda * M;
3099  M = q_mu * M;
3100  }
3101 
3102  // stress = 2*M*e(u) + L*tr(e(u))*I, where
3103  // e(u) = (1/2)*(grad(u) + grad(u)^T)
3104  const double M2 = 2.0*M;
3105  if (dim == 2)
3106  {
3107  L *= (grad(0,0) + grad(1,1));
3108  // order of the stress entries: s_xx, s_yy, s_xy
3109  flux(i+fnd*0) = M2*grad(0,0) + L;
3110  flux(i+fnd*1) = M2*grad(1,1) + L;
3111  flux(i+fnd*2) = M*(grad(0,1) + grad(1,0));
3112  }
3113  else if (dim == 3)
3114  {
3115  L *= (grad(0,0) + grad(1,1) + grad(2,2));
3116  // order of the stress entries: s_xx, s_yy, s_zz, s_xy, s_xz, s_yz
3117  flux(i+fnd*0) = M2*grad(0,0) + L;
3118  flux(i+fnd*1) = M2*grad(1,1) + L;
3119  flux(i+fnd*2) = M2*grad(2,2) + L;
3120  flux(i+fnd*3) = M*(grad(0,1) + grad(1,0));
3121  flux(i+fnd*4) = M*(grad(0,2) + grad(2,0));
3122  flux(i+fnd*5) = M*(grad(1,2) + grad(2,1));
3123  }
3124  }
3125 }
3126 
3127 double ElasticityIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
3129  Vector &flux, Vector *d_energy)
3130 {
3131  const int dof = fluxelem.GetDof();
3132  const int dim = fluxelem.GetDim();
3133  const int tdim = dim*(dim+1)/2; // num. entries in a symmetric tensor
3134  double L, M;
3135 
3136  // The MFEM_ASSERT constraints in ElasticityIntegrator::ComputeElementFlux
3137  // are assumed here too.
3138  MFEM_ASSERT(d_energy == NULL, "anisotropic estimates are not supported");
3139  MFEM_ASSERT(flux.Size() == dof*tdim, "invalid 'flux' vector");
3140 
3141 #ifndef MFEM_THREAD_SAFE
3142  shape.SetSize(dof);
3143 #else
3144  Vector shape(dof);
3145 #endif
3146  double pointstress_data[6];
3147  Vector pointstress(pointstress_data, tdim);
3148 
3149  // View of the 'flux' vector as a (dof x tdim) matrix
3150  DenseMatrix flux_mat(flux.GetData(), dof, tdim);
3151 
3152  // Use the same integration rule as in AssembleElementMatrix, replacing 'el'
3153  // with 'fluxelem' when 'IntRule' is not set.
3154  // Should we be using a different (more accurate) rule here?
3155  const IntegrationRule *ir = IntRule;
3156  if (ir == NULL)
3157  {
3158  int order = 2 * Trans.OrderGrad(&fluxelem);
3159  ir = &IntRules.Get(fluxelem.GetGeomType(), order);
3160  }
3161 
3162  double energy = 0.0;
3163 
3164  for (int i = 0; i < ir->GetNPoints(); i++)
3165  {
3166  const IntegrationPoint &ip = ir->IntPoint(i);
3167  fluxelem.CalcShape(ip, shape);
3168 
3169  flux_mat.MultTranspose(shape, pointstress);
3170 
3171  Trans.SetIntPoint(&ip);
3172  double w = Trans.Weight() * ip.weight;
3173 
3174  M = mu->Eval(Trans, ip);
3175  if (lambda)
3176  {
3177  L = lambda->Eval(Trans, ip);
3178  }
3179  else
3180  {
3181  L = q_lambda * M;
3182  M = q_mu * M;
3183  }
3184 
3185  // The strain energy density at a point is given by (1/2)*(s : e) where s
3186  // and e are the stress and strain tensors, respectively. Since we only
3187  // have the stress, we need to compute the strain from the stress:
3188  // s = 2*mu*e + lambda*tr(e)*I
3189  // Taking trace on both sides we find:
3190  // tr(s) = 2*mu*tr(e) + lambda*tr(e)*dim = (2*mu + dim*lambda)*tr(e)
3191  // which gives:
3192  // tr(e) = tr(s)/(2*mu + dim*lambda)
3193  // Then from the first identity above we can find the strain:
3194  // e = (1/(2*mu))*(s - lambda*tr(e)*I)
3195 
3196  double pt_e; // point strain energy density
3197  const double *s = pointstress_data;
3198  if (dim == 2)
3199  {
3200  // s entries: s_xx, s_yy, s_xy
3201  const double tr_e = (s[0] + s[1])/(2*(M + L));
3202  L *= tr_e;
3203  pt_e = (0.25/M)*(s[0]*(s[0] - L) + s[1]*(s[1] - L) + 2*s[2]*s[2]);
3204  }
3205  else // (dim == 3)
3206  {
3207  // s entries: s_xx, s_yy, s_zz, s_xy, s_xz, s_yz
3208  const double tr_e = (s[0] + s[1] + s[2])/(2*M + 3*L);
3209  L *= tr_e;
3210  pt_e = (0.25/M)*(s[0]*(s[0] - L) + s[1]*(s[1] - L) + s[2]*(s[2] - L) +
3211  2*(s[3]*s[3] + s[4]*s[4] + s[5]*s[5]));
3212  }
3213 
3214  energy += w * pt_e;
3215  }
3216 
3217  return energy;
3218 }
3219 
3220 void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
3221  const FiniteElement &el2,
3223  DenseMatrix &elmat)
3224 {
3225  int ndof1, ndof2;
3226 
3227  double un, a, b, w;
3228 
3229  dim = el1.GetDim();
3230  ndof1 = el1.GetDof();
3231  Vector vu(dim), nor(dim);
3232 
3233  if (Trans.Elem2No >= 0)
3234  {
3235  ndof2 = el2.GetDof();
3236  }
3237  else
3238  {
3239  ndof2 = 0;
3240  }
3241 
3242  shape1.SetSize(ndof1);
3243  shape2.SetSize(ndof2);
3244  elmat.SetSize(ndof1 + ndof2);
3245  elmat = 0.0;
3246 
3247  const IntegrationRule *ir = IntRule;
3248  if (ir == NULL)
3249  {
3250  int order;
3251  // Assuming order(u)==order(mesh)
3252  if (Trans.Elem2No >= 0)
3253  order = (min(Trans.Elem1->OrderW(), Trans.Elem2->OrderW()) +
3254  2*max(el1.GetOrder(), el2.GetOrder()));
3255  else
3256  {
3257  order = Trans.Elem1->OrderW() + 2*el1.GetOrder();
3258  }
3259  if (el1.Space() == FunctionSpace::Pk)
3260  {
3261  order++;
3262  }
3263  ir = &IntRules.Get(Trans.GetGeometryType(), order);
3264  }
3265 
3266  for (int p = 0; p < ir->GetNPoints(); p++)
3267  {
3268  const IntegrationPoint &ip = ir->IntPoint(p);
3269 
3270  // Set the integration point in the face and the neighboring elements
3271  Trans.SetAllIntPoints(&ip);
3272 
3273  // Access the neighboring elements' integration points
3274  // Note: eip2 will only contain valid data if Elem2 exists
3275  const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
3276  const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
3277 
3278  el1.CalcShape(eip1, shape1);
3279 
3280  u->Eval(vu, *Trans.Elem1, eip1);
3281 
3282  if (dim == 1)
3283  {
3284  nor(0) = 2*eip1.x - 1.0;
3285  }
3286  else
3287  {
3288  CalcOrtho(Trans.Jacobian(), nor);
3289  }
3290 
3291  un = vu * nor;
3292  a = 0.5 * alpha * un;
3293  b = beta * fabs(un);
3294  // note: if |alpha/2|==|beta| then |a|==|b|, i.e. (a==b) or (a==-b)
3295  // and therefore two blocks in the element matrix contribution
3296  // (from the current quadrature point) are 0
3297 
3298  if (rho)
3299  {
3300  double rho_p;
3301  if (un >= 0.0 && ndof2)
3302  {
3303  rho_p = rho->Eval(*Trans.Elem2, eip2);
3304  }
3305  else
3306  {
3307  rho_p = rho->Eval(*Trans.Elem1, eip1);
3308  }
3309  a *= rho_p;
3310  b *= rho_p;
3311  }
3312 
3313  w = ip.weight * (a+b);
3314  if (w != 0.0)
3315  {
3316  for (int i = 0; i < ndof1; i++)
3317  for (int j = 0; j < ndof1; j++)
3318  {
3319  elmat(i, j) += w * shape1(i) * shape1(j);
3320  }
3321  }
3322 
3323  if (ndof2)
3324  {
3325  el2.CalcShape(eip2, shape2);
3326 
3327  if (w != 0.0)
3328  for (int i = 0; i < ndof2; i++)
3329  for (int j = 0; j < ndof1; j++)
3330  {
3331  elmat(ndof1+i, j) -= w * shape2(i) * shape1(j);
3332  }
3333 
3334  w = ip.weight * (b-a);
3335  if (w != 0.0)
3336  {
3337  for (int i = 0; i < ndof2; i++)
3338  for (int j = 0; j < ndof2; j++)
3339  {
3340  elmat(ndof1+i, ndof1+j) += w * shape2(i) * shape2(j);
3341  }
3342 
3343  for (int i = 0; i < ndof1; i++)
3344  for (int j = 0; j < ndof2; j++)
3345  {
3346  elmat(i, ndof1+j) -= w * shape1(i) * shape2(j);
3347  }
3348  }
3349  }
3350  }
3351 }
3352 
3353 
3355  Geometry::Type geom, int order, FaceElementTransformations &T)
3356 {
3357  int int_order = T.Elem1->OrderW() + 2*order;
3358  return IntRules.Get(geom, int_order);
3359 }
3360 
3361 void DGDiffusionIntegrator::AssembleFaceMatrix(
3362  const FiniteElement &el1, const FiniteElement &el2,
3364 {
3365  int dim, ndof1, ndof2, ndofs;
3366  bool kappa_is_nonzero = (kappa != 0.);
3367  double w, wq = 0.0;
3368 
3369  dim = el1.GetDim();
3370  ndof1 = el1.GetDof();
3371 
3372  nor.SetSize(dim);
3373  nh.SetSize(dim);
3374  ni.SetSize(dim);
3375  adjJ.SetSize(dim);
3376  if (MQ)
3377  {
3378  mq.SetSize(dim);
3379  }
3380 
3381  shape1.SetSize(ndof1);
3382  dshape1.SetSize(ndof1, dim);
3383  dshape1dn.SetSize(ndof1);
3384  if (Trans.Elem2No >= 0)
3385  {
3386  ndof2 = el2.GetDof();
3387  shape2.SetSize(ndof2);
3388  dshape2.SetSize(ndof2, dim);
3389  dshape2dn.SetSize(ndof2);
3390  }
3391  else
3392  {
3393  ndof2 = 0;
3394  }
3395 
3396  ndofs = ndof1 + ndof2;
3397  elmat.SetSize(ndofs);
3398  elmat = 0.0;
3399  if (kappa_is_nonzero)
3400  {
3401  jmat.SetSize(ndofs);
3402  jmat = 0.;
3403  }
3404 
3405  const IntegrationRule *ir = IntRule;
3406  if (ir == NULL)
3407  {
3408  // a simple choice for the integration order; is this OK?
3409  int order;
3410  if (ndof2)
3411  {
3412  order = 2*max(el1.GetOrder(), el2.GetOrder());
3413  }
3414  else
3415  {
3416  order = 2*el1.GetOrder();
3417  }
3418  ir = &IntRules.Get(Trans.GetGeometryType(), order);
3419  }
3420 
3421  // assemble: < {(Q \nabla u).n},[v] > --> elmat
3422  // kappa < {h^{-1} Q} [u],[v] > --> jmat
3423  for (int p = 0; p < ir->GetNPoints(); p++)
3424  {
3425  const IntegrationPoint &ip = ir->IntPoint(p);
3426 
3427  // Set the integration point in the face and the neighboring elements
3428  Trans.SetAllIntPoints(&ip);
3429 
3430  // Access the neighboring elements' integration points
3431  // Note: eip2 will only contain valid data if Elem2 exists
3432  const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
3433  const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
3434 
3435  if (dim == 1)
3436  {
3437  nor(0) = 2*eip1.x - 1.0;
3438  }
3439  else
3440  {
3441  CalcOrtho(Trans.Jacobian(), nor);
3442  }
3443 
3444  el1.CalcShape(eip1, shape1);
3445  el1.CalcDShape(eip1, dshape1);
3446  w = ip.weight/Trans.Elem1->Weight();
3447  if (ndof2)
3448  {
3449  w /= 2;
3450  }
3451  if (!MQ)
3452  {
3453  if (Q)
3454  {
3455  w *= Q->Eval(*Trans.Elem1, eip1);
3456  }
3457  ni.Set(w, nor);
3458  }
3459  else
3460  {
3461  nh.Set(w, nor);
3462  MQ->Eval(mq, *Trans.Elem1, eip1);
3463  mq.MultTranspose(nh, ni);
3464  }
3465  CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
3466  adjJ.Mult(ni, nh);
3467  if (kappa_is_nonzero)
3468  {
3469  wq = ni * nor;
3470  }
3471  // Note: in the jump term, we use 1/h1 = |nor|/det(J1) which is
3472  // independent of Loc1 and always gives the size of element 1 in
3473  // direction perpendicular to the face. Indeed, for linear transformation
3474  // |nor|=measure(face)/measure(ref. face),
3475  // det(J1)=measure(element)/measure(ref. element),
3476  // and the ratios measure(ref. element)/measure(ref. face) are
3477  // compatible for all element/face pairs.
3478  // For example: meas(ref. tetrahedron)/meas(ref. triangle) = 1/3, and
3479  // for any tetrahedron vol(tet)=(1/3)*height*area(base).
3480  // For interior faces: q_e/h_e=(q1/h1+q2/h2)/2.
3481 
3482  dshape1.Mult(nh, dshape1dn);
3483  for (int i = 0; i < ndof1; i++)
3484  for (int j = 0; j < ndof1; j++)
3485  {
3486  elmat(i, j) += shape1(i) * dshape1dn(j);
3487  }
3488 
3489  if (ndof2)
3490  {
3491  el2.CalcShape(eip2, shape2);
3492  el2.CalcDShape(eip2, dshape2);
3493  w = ip.weight/2/Trans.Elem2->Weight();
3494  if (!MQ)
3495  {
3496  if (Q)
3497  {
3498  w *= Q->Eval(*Trans.Elem2, eip2);
3499  }
3500  ni.Set(w, nor);
3501  }
3502  else
3503  {
3504  nh.Set(w, nor);
3505  MQ->Eval(mq, *Trans.Elem2, eip2);
3506  mq.MultTranspose(nh, ni);
3507  }
3508  CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
3509  adjJ.Mult(ni, nh);
3510  if (kappa_is_nonzero)
3511  {
3512  wq += ni * nor;
3513  }
3514 
3515  dshape2.Mult(nh, dshape2dn);
3516 
3517  for (int i = 0; i < ndof1; i++)
3518  for (int j = 0; j < ndof2; j++)
3519  {
3520  elmat(i, ndof1 + j) += shape1(i) * dshape2dn(j);
3521  }
3522 
3523  for (int i = 0; i < ndof2; i++)
3524  for (int j = 0; j < ndof1; j++)
3525  {
3526  elmat(ndof1 + i, j) -= shape2(i) * dshape1dn(j);
3527  }
3528 
3529  for (int i = 0; i < ndof2; i++)
3530  for (int j = 0; j < ndof2; j++)
3531  {
3532  elmat(ndof1 + i, ndof1 + j) -= shape2(i) * dshape2dn(j);
3533  }
3534  }
3535 
3536  if (kappa_is_nonzero)
3537  {
3538  // only assemble the lower triangular part of jmat
3539  wq *= kappa;
3540  for (int i = 0; i < ndof1; i++)
3541  {
3542  const double wsi = wq*shape1(i);
3543  for (int j = 0; j <= i; j++)
3544  {
3545  jmat(i, j) += wsi * shape1(j);
3546  }
3547  }
3548  if (ndof2)
3549  {
3550  for (int i = 0; i < ndof2; i++)
3551  {
3552  const int i2 = ndof1 + i;
3553  const double wsi = wq*shape2(i);
3554  for (int j = 0; j < ndof1; j++)
3555  {
3556  jmat(i2, j) -= wsi * shape1(j);
3557  }
3558  for (int j = 0; j <= i; j++)
3559  {
3560  jmat(i2, ndof1 + j) += wsi * shape2(j);
3561  }
3562  }
3563  }
3564  }
3565  }
3566 
3567  // elmat := -elmat + sigma*elmat^t + jmat
3568  if (kappa_is_nonzero)
3569  {
3570  for (int i = 0; i < ndofs; i++)
3571  {
3572  for (int j = 0; j < i; j++)
3573  {
3574  double aij = elmat(i,j), aji = elmat(j,i), mij = jmat(i,j);
3575  elmat(i,j) = sigma*aji - aij + mij;
3576  elmat(j,i) = sigma*aij - aji + mij;
3577  }
3578  elmat(i,i) = (sigma - 1.)*elmat(i,i) + jmat(i,i);
3579  }
3580  }
3581  else
3582  {
3583  for (int i = 0; i < ndofs; i++)
3584  {
3585  for (int j = 0; j < i; j++)
3586  {
3587  double aij = elmat(i,j), aji = elmat(j,i);
3588  elmat(i,j) = sigma*aji - aij;
3589  elmat(j,i) = sigma*aij - aji;
3590  }
3591  elmat(i,i) *= (sigma - 1.);
3592  }
3593  }
3594 }
3595 
3596 
3597 // static method
3598 void DGElasticityIntegrator::AssembleBlock(
3599  const int dim, const int row_ndofs, const int col_ndofs,
3600  const int row_offset, const int col_offset,
3601  const double jmatcoef, const Vector &col_nL, const Vector &col_nM,
3602  const Vector &row_shape, const Vector &col_shape,
3603  const Vector &col_dshape_dnM, const DenseMatrix &col_dshape,
3604  DenseMatrix &elmat, DenseMatrix &jmat)
3605 {
3606  for (int jm = 0, j = col_offset; jm < dim; ++jm)
3607  {
3608  for (int jdof = 0; jdof < col_ndofs; ++jdof, ++j)
3609  {
3610  const double t2 = col_dshape_dnM(jdof);
3611  for (int im = 0, i = row_offset; im < dim; ++im)
3612  {
3613  const double t1 = col_dshape(jdof, jm) * col_nL(im);
3614  const double t3 = col_dshape(jdof, im) * col_nM(jm);
3615  const double tt = t1 + ((im == jm) ? t2 : 0.0) + t3;
3616  for (int idof = 0; idof < row_ndofs; ++idof, ++i)
3617  {
3618  elmat(i, j) += row_shape(idof) * tt;
3619  }
3620  }
3621  }
3622  }
3623 
3624  if (jmatcoef == 0.0) { return; }
3625 
3626  for (int d = 0; d < dim; ++d)
3627  {
3628  const int jo = col_offset + d*col_ndofs;
3629  const int io = row_offset + d*row_ndofs;
3630  for (int jdof = 0, j = jo; jdof < col_ndofs; ++jdof, ++j)
3631  {
3632  const double sj = jmatcoef * col_shape(jdof);
3633  for (int i = max(io,j), idof = i - io; idof < row_ndofs; ++idof, ++i)
3634  {
3635  jmat(i, j) += row_shape(idof) * sj;
3636  }
3637  }
3638  }
3639 }
3640 
3641 void DGElasticityIntegrator::AssembleFaceMatrix(
3642  const FiniteElement &el1, const FiniteElement &el2,
3644 {
3645 #ifdef MFEM_THREAD_SAFE
3646  // For descriptions of these variables, see the class declaration.
3647  Vector shape1, shape2;
3648  DenseMatrix dshape1, dshape2;
3649  DenseMatrix adjJ;
3650  DenseMatrix dshape1_ps, dshape2_ps;
3651  Vector nor;
3652  Vector nL1, nL2;
3653  Vector nM1, nM2;
3654  Vector dshape1_dnM, dshape2_dnM;
3655  DenseMatrix jmat;
3656 #endif
3657 
3658  const int dim = el1.GetDim();
3659  const int ndofs1 = el1.GetDof();
3660  const int ndofs2 = (Trans.Elem2No >= 0) ? el2.GetDof() : 0;
3661  const int nvdofs = dim*(ndofs1 + ndofs2);
3662 
3663  // Initially 'elmat' corresponds to the term:
3664  // < { sigma(u) . n }, [v] > =
3665  // < { (lambda div(u) I + mu (grad(u) + grad(u)^T)) . n }, [v] >
3666  // But eventually, it's going to be replaced by:
3667  // elmat := -elmat + alpha*elmat^T + jmat
3668  elmat.SetSize(nvdofs);
3669  elmat = 0.;
3670 
3671  const bool kappa_is_nonzero = (kappa != 0.0);
3672  if (kappa_is_nonzero)
3673  {
3674  jmat.SetSize(nvdofs);
3675  jmat = 0.;
3676  }
3677 
3678  adjJ.SetSize(dim);
3679  shape1.SetSize(ndofs1);
3680  dshape1.SetSize(ndofs1, dim);
3681  dshape1_ps.SetSize(ndofs1, dim);
3682  nor.SetSize(dim);
3683  nL1.SetSize(dim);
3684  nM1.SetSize(dim);
3685  dshape1_dnM.SetSize(ndofs1);
3686 
3687  if (ndofs2)
3688  {
3689  shape2.SetSize(ndofs2);
3690  dshape2.SetSize(ndofs2, dim);
3691  dshape2_ps.SetSize(ndofs2, dim);
3692  nL2.SetSize(dim);
3693  nM2.SetSize(dim);
3694  dshape2_dnM.SetSize(ndofs2);
3695  }
3696 
3697  const IntegrationRule *ir = IntRule;
3698  if (ir == NULL)
3699  {
3700  // a simple choice for the integration order; is this OK?
3701  const int order = 2 * max(el1.GetOrder(), ndofs2 ? el2.GetOrder() : 0);
3702  ir = &IntRules.Get(Trans.GetGeometryType(), order);
3703  }
3704 
3705  for (int pind = 0; pind < ir->GetNPoints(); ++pind)
3706  {
3707  const IntegrationPoint &ip = ir->IntPoint(pind);
3708 
3709  // Set the integration point in the face and the neighboring elements
3710  Trans.SetAllIntPoints(&ip);
3711 
3712  // Access the neighboring elements' integration points
3713  // Note: eip2 will only contain valid data if Elem2 exists
3714  const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
3715  const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
3716 
3717  el1.CalcShape(eip1, shape1);
3718  el1.CalcDShape(eip1, dshape1);
3719 
3720  CalcAdjugate(Trans.Elem1->Jacobian(), adjJ);
3721  Mult(dshape1, adjJ, dshape1_ps);
3722 
3723  if (dim == 1)
3724  {
3725  nor(0) = 2*eip1.x - 1.0;
3726  }
3727  else
3728  {
3729  CalcOrtho(Trans.Jacobian(), nor);
3730  }
3731 
3732  double w, wLM;
3733  if (ndofs2)
3734  {
3735  el2.CalcShape(eip2, shape2);
3736  el2.CalcDShape(eip2, dshape2);
3737  CalcAdjugate(Trans.Elem2->Jacobian(), adjJ);
3738  Mult(dshape2, adjJ, dshape2_ps);
3739 
3740  w = ip.weight/2;
3741  const double w2 = w / Trans.Elem2->Weight();
3742  const double wL2 = w2 * lambda->Eval(*Trans.Elem2, eip2);
3743  const double wM2 = w2 * mu->Eval(*Trans.Elem2, eip2);
3744  nL2.Set(wL2, nor);
3745  nM2.Set(wM2, nor);
3746  wLM = (wL2 + 2.0*wM2);
3747  dshape2_ps.Mult(nM2, dshape2_dnM);
3748  }
3749  else
3750  {
3751  w = ip.weight;
3752  wLM = 0.0;
3753  }
3754 
3755  {
3756  const double w1 = w / Trans.Elem1->Weight();
3757  const double wL1 = w1 * lambda->Eval(*Trans.Elem1, eip1);
3758  const double wM1 = w1 * mu->Eval(*Trans.Elem1, eip1);
3759  nL1.Set(wL1, nor);
3760  nM1.Set(wM1, nor);
3761  wLM += (wL1 + 2.0*wM1);
3762  dshape1_ps.Mult(nM1, dshape1_dnM);
3763  }
3764 
3765  const double jmatcoef = kappa * (nor*nor) * wLM;
3766 
3767  // (1,1) block
3768  AssembleBlock(
3769  dim, ndofs1, ndofs1, 0, 0, jmatcoef, nL1, nM1,
3770  shape1, shape1, dshape1_dnM, dshape1_ps, elmat, jmat);
3771 
3772  if (ndofs2 == 0) { continue; }
3773 
3774  // In both elmat and jmat, shape2 appears only with a minus sign.
3775  shape2.Neg();
3776 
3777  // (1,2) block
3778  AssembleBlock(
3779  dim, ndofs1, ndofs2, 0, dim*ndofs1, jmatcoef, nL2, nM2,
3780  shape1, shape2, dshape2_dnM, dshape2_ps, elmat, jmat);
3781  // (2,1) block
3782  AssembleBlock(
3783  dim, ndofs2, ndofs1, dim*ndofs1, 0, jmatcoef, nL1, nM1,
3784  shape2, shape1, dshape1_dnM, dshape1_ps, elmat, jmat);
3785  // (2,2) block
3786  AssembleBlock(
3787  dim, ndofs2, ndofs2, dim*ndofs1, dim*ndofs1, jmatcoef, nL2, nM2,
3788  shape2, shape2, dshape2_dnM, dshape2_ps, elmat, jmat);
3789  }
3790 
3791  // elmat := -elmat + alpha*elmat^t + jmat
3792  if (kappa_is_nonzero)
3793  {
3794  for (int i = 0; i < nvdofs; ++i)
3795  {
3796  for (int j = 0; j < i; ++j)
3797  {
3798  double aij = elmat(i,j), aji = elmat(j,i), mij = jmat(i,j);
3799  elmat(i,j) = alpha*aji - aij + mij;
3800  elmat(j,i) = alpha*aij - aji + mij;
3801  }
3802  elmat(i,i) = (alpha - 1.)*elmat(i,i) + jmat(i,i);
3803  }
3804  }
3805  else
3806  {
3807  for (int i = 0; i < nvdofs; ++i)
3808  {
3809  for (int j = 0; j < i; ++j)
3810  {
3811  double aij = elmat(i,j), aji = elmat(j,i);
3812  elmat(i,j) = alpha*aji - aij;
3813  elmat(j,i) = alpha*aij - aji;
3814  }
3815  elmat(i,i) *= (alpha - 1.);
3816  }
3817  }
3818 }
3819 
3820 
3821 void TraceJumpIntegrator::AssembleFaceMatrix(
3822  const FiniteElement &trial_face_fe, const FiniteElement &test_fe1,
3824  DenseMatrix &elmat)
3825 {
3826  int i, j, face_ndof, ndof1, ndof2;
3827  int order;
3828 
3829  double w;
3830 
3831  face_ndof = trial_face_fe.GetDof();
3832  ndof1 = test_fe1.GetDof();
3833 
3834  face_shape.SetSize(face_ndof);
3835  shape1.SetSize(ndof1);
3836 
3837  if (Trans.Elem2No >= 0)
3838  {
3839  ndof2 = test_fe2.GetDof();
3840  shape2.SetSize(ndof2);
3841  }
3842  else
3843  {
3844  ndof2 = 0;
3845  }
3846 
3847  elmat.SetSize(ndof1 + ndof2, face_ndof);
3848  elmat = 0.0;
3849 
3850  const IntegrationRule *ir = IntRule;
3851  if (ir == NULL)
3852  {
3853  if (Trans.Elem2No >= 0)
3854  {
3855  order = max(test_fe1.GetOrder(), test_fe2.GetOrder());
3856  }
3857  else
3858  {
3859  order = test_fe1.GetOrder();
3860  }
3861  order += trial_face_fe.GetOrder();
3862  if (trial_face_fe.GetMapType() == FiniteElement::VALUE)
3863  {
3864  order += Trans.OrderW();
3865  }
3866  ir = &IntRules.Get(Trans.GetGeometryType(), order);
3867  }
3868 
3869  for (int p = 0; p < ir->GetNPoints(); p++)
3870  {
3871  const IntegrationPoint &ip = ir->IntPoint(p);
3872 
3873  // Set the integration point in the face and the neighboring elements
3874  Trans.SetAllIntPoints(&ip);
3875 
3876  // Access the neighboring elements' integration points
3877  // Note: eip2 will only contain valid data if Elem2 exists
3878  const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
3879  const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
3880 
3881  // Trace finite element shape function
3882  trial_face_fe.CalcShape(ip, face_shape);
3883  // Side 1 finite element shape function
3884  test_fe1.CalcShape(eip1, shape1);
3885  if (ndof2)
3886  {
3887  // Side 2 finite element shape function
3888  test_fe2.CalcShape(eip2, shape2);
3889  }
3890  w = ip.weight;
3891  if (trial_face_fe.GetMapType() == FiniteElement::VALUE)
3892  {
3893  w *= Trans.Weight();
3894  }
3895  face_shape *= w;
3896  for (i = 0; i < ndof1; i++)
3897  for (j = 0; j < face_ndof; j++)
3898  {
3899  elmat(i, j) += shape1(i) * face_shape(j);
3900  }
3901  if (ndof2)
3902  {
3903  // Subtract contribution from side 2
3904  for (i = 0; i < ndof2; i++)
3905  for (j = 0; j < face_ndof; j++)
3906  {
3907  elmat(ndof1+i, j) -= shape2(i) * face_shape(j);
3908  }
3909  }
3910  }
3911 }
3912 
3913 void NormalTraceJumpIntegrator::AssembleFaceMatrix(
3914  const FiniteElement &trial_face_fe, const FiniteElement &test_fe1,
3916  DenseMatrix &elmat)
3917 {
3918  int i, j, face_ndof, ndof1, ndof2, dim;
3919  int order;
3920 
3921  MFEM_VERIFY(trial_face_fe.GetMapType() == FiniteElement::VALUE, "");
3922 
3923  face_ndof = trial_face_fe.GetDof();
3924  ndof1 = test_fe1.GetDof();
3925  dim = test_fe1.GetDim();
3926 
3927  face_shape.SetSize(face_ndof);
3928  normal.SetSize(dim);
3929  shape1.SetSize(ndof1,dim);
3930  shape1_n.SetSize(ndof1);
3931 
3932  if (Trans.Elem2No >= 0)
3933  {
3934  ndof2 = test_fe2.GetDof();
3935  shape2.SetSize(ndof2,dim);
3936  shape2_n.SetSize(ndof2);
3937  }
3938  else
3939  {
3940  ndof2 = 0;
3941  }
3942 
3943  elmat.SetSize(ndof1 + ndof2, face_ndof);
3944  elmat = 0.0;
3945 
3946  const IntegrationRule *ir = IntRule;
3947  if (ir == NULL)
3948  {
3949  if (Trans.Elem2No >= 0)
3950  {
3951  order = max(test_fe1.GetOrder(), test_fe2.GetOrder()) - 1;
3952  }
3953  else
3954  {
3955  order = test_fe1.GetOrder() - 1;
3956  }
3957  order += trial_face_fe.GetOrder();
3958  ir = &IntRules.Get(Trans.GetGeometryType(), order);
3959  }
3960 
3961  for (int p = 0; p < ir->GetNPoints(); p++)
3962  {
3963  const IntegrationPoint &ip = ir->IntPoint(p);
3964  IntegrationPoint eip1, eip2;
3965  // Trace finite element shape function
3966  trial_face_fe.CalcShape(ip, face_shape);
3967  Trans.Loc1.Transf.SetIntPoint(&ip);
3968  CalcOrtho(Trans.Loc1.Transf.Jacobian(), normal);
3969  // Side 1 finite element shape function
3970  Trans.Loc1.Transform(ip, eip1);
3971  test_fe1.CalcVShape(eip1, shape1);
3972  shape1.Mult(normal, shape1_n);
3973  if (ndof2)
3974  {
3975  // Side 2 finite element shape function
3976  Trans.Loc2.Transform(ip, eip2);
3977  test_fe2.CalcVShape(eip2, shape2);
3978  Trans.Loc2.Transf.SetIntPoint(&ip);
3979  CalcOrtho(Trans.Loc2.Transf.Jacobian(), normal);
3980  shape2.Mult(normal, shape2_n);
3981  }
3982  face_shape *= ip.weight;
3983  for (i = 0; i < ndof1; i++)
3984  for (j = 0; j < face_ndof; j++)
3985  {
3986  elmat(i, j) += shape1_n(i) * face_shape(j);
3987  }
3988  if (ndof2)
3989  {
3990  // Subtract contribution from side 2
3991  for (i = 0; i < ndof2; i++)
3992  for (j = 0; j < face_ndof; j++)
3993  {
3994  elmat(ndof1+i, j) -= shape2_n(i) * face_shape(j);
3995  }
3996  }
3997  }
3998 }
3999 
4000 
4001 void NormalInterpolator::AssembleElementMatrix2(
4002  const FiniteElement &dom_fe, const FiniteElement &ran_fe,
4004 {
4005  int spaceDim = Trans.GetSpaceDim();
4006  elmat.SetSize(ran_fe.GetDof(), spaceDim*dom_fe.GetDof());
4007  Vector n(spaceDim), shape(dom_fe.GetDof());
4008 
4009  const IntegrationRule &ran_nodes = ran_fe.GetNodes();
4010  for (int i = 0; i < ran_nodes.Size(); i++)
4011  {
4012  const IntegrationPoint &ip = ran_nodes.IntPoint(i);
4013  Trans.SetIntPoint(&ip);
4014  CalcOrtho(Trans.Jacobian(), n);
4015  dom_fe.CalcShape(ip, shape);
4016  for (int j = 0; j < shape.Size(); j++)
4017  {
4018  for (int d = 0; d < spaceDim; d++)
4019  {
4020  elmat(i, j+d*shape.Size()) = shape(j)*n(d);
4021  }
4022  }
4023  }
4024 }
4025 
4026 
4027 namespace internal
4028 {
4029 
4030 // Scalar shape functions scaled by scalar coefficient.
4031 // Used in the implementation of class ScalarProductInterpolator below.
4032 struct ShapeCoefficient : public VectorCoefficient
4033 {
4034  Coefficient &Q;
4035  const FiniteElement &fe;
4036 
4037  ShapeCoefficient(Coefficient &q, const FiniteElement &fe_)
4038  : VectorCoefficient(fe_.GetDof()), Q(q), fe(fe_) { }
4039 
4040  using VectorCoefficient::Eval;
4041  virtual void Eval(Vector &V, ElementTransformation &T,
4042  const IntegrationPoint &ip)
4043  {
4044  V.SetSize(vdim);
4045  fe.CalcPhysShape(T, V);
4046  V *= Q.Eval(T, ip);
4047  }
4048 };
4049 
4050 }
4051 
4052 void
4053 ScalarProductInterpolator::AssembleElementMatrix2(const FiniteElement &dom_fe,
4054  const FiniteElement &ran_fe,
4056  DenseMatrix &elmat)
4057 {
4058  internal::ShapeCoefficient dom_shape_coeff(*Q, dom_fe);
4059 
4060  elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
4061 
4062  Vector elmat_as_vec(elmat.Data(), ran_fe.GetDof()*dom_fe.GetDof());
4063 
4064  ran_fe.Project(dom_shape_coeff, Trans, elmat_as_vec);
4065 }
4066 
4067 
4068 void
4069 ScalarVectorProductInterpolator::AssembleElementMatrix2(
4070  const FiniteElement &dom_fe,
4071  const FiniteElement &ran_fe,
4073  DenseMatrix &elmat)
4074 {
4075  // Vector shape functions scaled by scalar coefficient
4076  struct VShapeCoefficient : public MatrixCoefficient
4077  {
4078  Coefficient &Q;
4079  const FiniteElement &fe;
4080 
4081  VShapeCoefficient(Coefficient &q, const FiniteElement &fe_, int sdim)
4082  : MatrixCoefficient(fe_.GetDof(), sdim), Q(q), fe(fe_) { }
4083 
4084  virtual void Eval(DenseMatrix &M, ElementTransformation &T,
4085  const IntegrationPoint &ip)
4086  {
4087  M.SetSize(height, width);
4088  fe.CalcPhysVShape(T, M);
4089  M *= Q.Eval(T, ip);
4090  }
4091  };
4092 
4093  VShapeCoefficient dom_shape_coeff(*Q, dom_fe, Trans.GetSpaceDim());
4094 
4095  elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
4096 
4097  Vector elmat_as_vec(elmat.Data(), ran_fe.GetDof()*dom_fe.GetDof());
4098 
4099  ran_fe.ProjectMatrixCoefficient(dom_shape_coeff, Trans, elmat_as_vec);
4100 }
4101 
4102 
4103 void
4104 VectorScalarProductInterpolator::AssembleElementMatrix2(
4105  const FiniteElement &dom_fe,
4106  const FiniteElement &ran_fe,
4108  DenseMatrix &elmat)
4109 {
4110  // Scalar shape functions scaled by vector coefficient
4111  struct VecShapeCoefficient : public MatrixCoefficient
4112  {
4113  VectorCoefficient &VQ;
4114  const FiniteElement &fe;
4115  Vector vc, shape;
4116 
4117  VecShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
4118  : MatrixCoefficient(fe_.GetDof(), vq.GetVDim()), VQ(vq), fe(fe_),
4119  vc(width), shape(height) { }
4120 
4121  virtual void Eval(DenseMatrix &M, ElementTransformation &T,
4122  const IntegrationPoint &ip)
4123  {
4124  M.SetSize(height, width);
4125  VQ.Eval(vc, T, ip);
4126  fe.CalcPhysShape(T, shape);
4127  MultVWt(shape, vc, M);
4128  }
4129  };
4130 
4131  VecShapeCoefficient dom_shape_coeff(*VQ, dom_fe);
4132 
4133  elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
4134 
4135  Vector elmat_as_vec(elmat.Data(), ran_fe.GetDof()*dom_fe.GetDof());
4136 
4137  ran_fe.ProjectMatrixCoefficient(dom_shape_coeff, Trans, elmat_as_vec);
4138 }
4139 
4140 
4141 void
4142 ScalarCrossProductInterpolator::AssembleElementMatrix2(
4143  const FiniteElement &dom_fe,
4144  const FiniteElement &ran_fe,
4146  DenseMatrix &elmat)
4147 {
4148  // Vector coefficient product with vector shape functions
4149  struct VCrossVShapeCoefficient : public VectorCoefficient
4150  {
4151  VectorCoefficient &VQ;
4152  const FiniteElement &fe;
4153  DenseMatrix vshape;
4154  Vector vc;
4155 
4156  VCrossVShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
4157  : VectorCoefficient(fe_.GetDof()), VQ(vq), fe(fe_),
4158  vshape(vdim, vq.GetVDim()), vc(vq.GetVDim()) { }
4159 
4160  using VectorCoefficient::Eval;
4161  virtual void Eval(Vector &V, ElementTransformation &T,
4162  const IntegrationPoint &ip)
4163  {
4164  V.SetSize(vdim);
4165  VQ.Eval(vc, T, ip);
4166  fe.CalcPhysVShape(T, vshape);
4167  for (int k = 0; k < vdim; k++)
4168  {
4169  V(k) = vc(0) * vshape(k,1) - vc(1) * vshape(k,0);
4170  }
4171  }
4172  };
4173 
4174  VCrossVShapeCoefficient dom_shape_coeff(*VQ, dom_fe);
4175 
4176  elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
4177 
4178  Vector elmat_as_vec(elmat.Data(), elmat.Height()*elmat.Width());
4179 
4180  ran_fe.Project(dom_shape_coeff, Trans, elmat_as_vec);
4181 }
4182 
4183 void
4184 VectorCrossProductInterpolator::AssembleElementMatrix2(
4185  const FiniteElement &dom_fe,
4186  const FiniteElement &ran_fe,
4188  DenseMatrix &elmat)
4189 {
4190  // Vector coefficient product with vector shape functions
4191  struct VCrossVShapeCoefficient : public MatrixCoefficient
4192  {
4193  VectorCoefficient &VQ;
4194  const FiniteElement &fe;
4195  DenseMatrix vshape;
4196  Vector vc;
4197 
4198  VCrossVShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
4199  : MatrixCoefficient(fe_.GetDof(), vq.GetVDim()), VQ(vq), fe(fe_),
4200  vshape(height, width), vc(width)
4201  {
4202  MFEM_ASSERT(width == 3, "");
4203  }
4204 
4205  virtual void Eval(DenseMatrix &M, ElementTransformation &T,
4206  const IntegrationPoint &ip)
4207  {
4208  M.SetSize(height, width);
4209  VQ.Eval(vc, T, ip);
4210  fe.CalcPhysVShape(T, vshape);
4211  for (int k = 0; k < height; k++)
4212  {
4213  M(k,0) = vc(1) * vshape(k,2) - vc(2) * vshape(k,1);
4214  M(k,1) = vc(2) * vshape(k,0) - vc(0) * vshape(k,2);
4215  M(k,2) = vc(0) * vshape(k,1) - vc(1) * vshape(k,0);
4216  }
4217  }
4218  };
4219 
4220  VCrossVShapeCoefficient dom_shape_coeff(*VQ, dom_fe);
4221 
4222  if (ran_fe.GetRangeType() == FiniteElement::SCALAR)
4223  {
4224  elmat.SetSize(ran_fe.GetDof()*VQ->GetVDim(),dom_fe.GetDof());
4225  }
4226  else
4227  {
4228  elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
4229  }
4230 
4231  Vector elmat_as_vec(elmat.Data(), elmat.Height()*elmat.Width());
4232 
4233  ran_fe.ProjectMatrixCoefficient(dom_shape_coeff, Trans, elmat_as_vec);
4234 }
4235 
4236 
4237 namespace internal
4238 {
4239 
4240 // Vector shape functions dot product with a vector coefficient.
4241 // Used in the implementation of class VectorInnerProductInterpolator below.
4242 struct VDotVShapeCoefficient : public VectorCoefficient
4243 {
4244  VectorCoefficient &VQ;
4245  const FiniteElement &fe;
4246  DenseMatrix vshape;
4247  Vector vc;
4248 
4249  VDotVShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
4250  : VectorCoefficient(fe_.GetDof()), VQ(vq), fe(fe_),
4251  vshape(vdim, vq.GetVDim()), vc(vq.GetVDim()) { }
4252 
4253  using VectorCoefficient::Eval;
4254  virtual void Eval(Vector &V, ElementTransformation &T,
4255  const IntegrationPoint &ip)
4256  {
4257  V.SetSize(vdim);
4258  VQ.Eval(vc, T, ip);
4259  fe.CalcPhysVShape(T, vshape);
4260  vshape.Mult(vc, V);
4261  }
4262 };
4263 
4264 }
4265 
4266 void
4267 VectorInnerProductInterpolator::AssembleElementMatrix2(
4268  const FiniteElement &dom_fe,
4269  const FiniteElement &ran_fe,
4271  DenseMatrix &elmat)
4272 {
4273  internal::VDotVShapeCoefficient dom_shape_coeff(*VQ, dom_fe);
4274 
4275  elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
4276 
4277  Vector elmat_as_vec(elmat.Data(), elmat.Height()*elmat.Width());
4278 
4279  ran_fe.Project(dom_shape_coeff, Trans, elmat_as_vec);
4280 }
4281 
4282 }
Abstract class for all finite elements.
Definition: fe_base.hpp:232
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:2781
const DenseMatrix & AdjugateJacobian()
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition: eltrans.hpp:135
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
Definition: densemat.cpp:3166
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual int OrderGrad(const FiniteElement *fe) const =0
Return the order of .
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:162
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
Definition: densemat.cpp:3147
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
Base class for vector Coefficients that optionally depend on time and space.
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe_base.cpp:64
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe_base.hpp:336
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
Definition: fe_base.hpp:316
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2479
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre)
A global object with all refined integration rules.
Definition: intrules.hpp:382
void CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape) const
Equivalent to the CalcVShape() method with the same arguments.
Definition: fe_base.hpp:410
void AddMult_a_ABt(double a, const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += a * A * B^t.
Definition: densemat.cpp:2979
double kappa
Definition: ex24.cpp:54
STL namespace.
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2693
void Mult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt = a * A * A^t.
Definition: densemat.cpp:3120
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:314
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_base.cpp:125
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
void Reset(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:90
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:215
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:154
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:3211
double b
Definition: lissajous.cpp:42
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
double InnerProduct(const double *x, const double *y) const
Compute y^t A x.
Definition: densemat.cpp:336
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...
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:679
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition: fe_base.cpp:70
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:3233
int GetVDim()
Returns dimension of the vector.
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:167
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition: fe_base.cpp:191
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2587
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Definition: fe_base.cpp:143
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1420
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:310
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
Definition: densemat.cpp:3136
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
Definition: densemat.cpp:2886
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void ClearExternalData()
Definition: densemat.hpp:95
Base class for Matrix Coefficients that optionally depend on time and space.
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2719
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0<=i<A.Height, 0<=j<A.Width.
Definition: densemat.cpp:1723
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
Vector & Set(const double a, const Vector &x)
(*this) = a * x
Definition: vector.cpp:268
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
const IntegrationRule & GetRule(const Integrator &integ, const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans)
double a
Definition: lissajous.cpp:41
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_base.cpp:181
void AddMultADBt(const DenseMatrix &A, const Vector &D, const DenseMatrix &B, DenseMatrix &ADBt)
ADBt = A D B^t, where D is diagonal.
Definition: densemat.cpp:2943
Class for integration point with weight.
Definition: intrules.hpp:25
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
Definition: fe_base.hpp:313
double mu
Definition: ex25.cpp:139
int dim
Definition: ex24.cpp:53
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_base.cpp:39
void GradToCurl(DenseMatrix &curl)
Definition: densemat.cpp:1478
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:348
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
const double alpha
Definition: ex15.cpp:369
Vector data type.
Definition: vector.hpp:60
void AddMultADAt(const DenseMatrix &A, const Vector &D, DenseMatrix &ADAt)
ADAt += A D A^t, where D is diagonal.
Definition: densemat.cpp:2737
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe_base.hpp:339
RefCoord s[3]
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:388
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...
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:326
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:3095
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Definition: fe_base.cpp:51
void Neg()
(*this) = -(*this)
Definition: vector.cpp:306
double sigma(const Vector &x)
Definition: maxwell.cpp:102