MFEM  v4.5.2
Finite element discretization library
coefficient.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 Coefficient class
13 
14 #include "fem.hpp"
15 
16 #include <cmath>
17 #include <limits>
18 
19 namespace mfem
20 {
21 
22 using namespace std;
23 
24 // Given an ElementTransformation and IntegrationPoint in a refined mesh,
25 // return the ElementTransformation of the parent coarse element, and set
26 // coarse_ip to the location of the original ip within the coarse element.
28  Mesh &coarse_mesh, const ElementTransformation &T,
29  const IntegrationPoint &ip, IntegrationPoint &coarse_ip)
30 {
31  Mesh &fine_mesh = *T.mesh;
32  // Get the element transformation of the coarse element containing the
33  // fine element.
34  int fine_element = T.ElementNo;
35  const CoarseFineTransformations &cf = fine_mesh.GetRefinementTransforms();
36  int coarse_element = cf.embeddings[fine_element].parent;
37  ElementTransformation *coarse_T = coarse_mesh.GetElementTransformation(
38  coarse_element);
39  // Transform the integration point from fine element coordinates to coarse
40  // element coordinates.
41  Geometry::Type geom = T.GetGeometryType();
42  IntegrationPointTransformation fine_to_coarse;
43  IsoparametricTransformation &emb_tr = fine_to_coarse.Transf;
44  emb_tr.SetIdentityTransformation(geom);
45  emb_tr.SetPointMat(cf.point_matrices[geom](cf.embeddings[fine_element].matrix));
46  fine_to_coarse.Transform(ip, coarse_ip);
47  coarse_T->SetIntPoint(&coarse_ip);
48  return coarse_T;
49 }
50 
52 {
53  QuadratureSpaceBase &qspace = *qf.GetSpace();
54  const int ne = qspace.GetNE();
55  Vector values;
56  for (int iel = 0; iel < ne; ++iel)
57  {
58  qf.GetValues(iel, values);
59  const IntegrationRule &ir = qspace.GetIntRule(iel);
60  ElementTransformation& T = *qspace.GetTransformation(iel);
61  for (int iq = 0; iq < ir.Size(); ++iq)
62  {
63  const IntegrationPoint &ip = ir[iq];
64  T.SetIntPoint(&ip);
65  const int iq_p = qspace.GetPermutedIndex(iel, iq);
66  values[iq_p] = Eval(T, ip);
67  }
68  }
69 }
70 
72 {
73  qf = constant;
74 }
75 
77  const IntegrationPoint & ip)
78 {
79  int att = T.Attribute;
80  return (constants(att-1));
81 }
82 
83 void PWCoefficient::InitMap(const Array<int> & attr,
84  const Array<Coefficient*> & coefs)
85 {
86  MFEM_VERIFY(attr.Size() == coefs.Size(),
87  "PWCoefficient: "
88  "Attribute and coefficient arrays have incompatible "
89  "dimensions.");
90 
91  for (int i=0; i<attr.Size(); i++)
92  {
93  if (coefs[i] != NULL)
94  {
95  UpdateCoefficient(attr[i], *coefs[i]);
96  }
97  }
98 }
99 
101 {
103 
104  std::map<int, Coefficient*>::iterator p = pieces.begin();
105  for (; p != pieces.end(); p++)
106  {
107  if (p->second != NULL)
108  {
109  p->second->SetTime(t);
110  }
111  }
112 }
113 
115  const IntegrationPoint &ip)
116 {
117  const int att = T.Attribute;
118  std::map<int, Coefficient*>::const_iterator p = pieces.find(att);
119  if (p != pieces.end())
120  {
121  if ( p->second != NULL)
122  {
123  return p->second->Eval(T, ip);
124  }
125  }
126  return 0.0;
127 }
128 
130  const IntegrationPoint & ip)
131 {
132  double x[3];
133  Vector transip(x, 3);
134 
135  T.Transform(ip, transip);
136 
137  if (Function)
138  {
139  return Function(transip);
140  }
141  else
142  {
143  return TDFunction(transip, GetTime());
144  }
145 }
146 
148  const IntegrationPoint &ip)
149 {
150  Mesh *gf_mesh = GridF->FESpace()->GetMesh();
151  if (T.mesh == gf_mesh)
152  {
153  return GridF->GetValue(T, ip, Component);
154  }
155  else
156  {
157  IntegrationPoint coarse_ip;
158  ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
159  return GridF->GetValue(*coarse_T, coarse_ip, Component);
160  }
161 }
162 
164 {
165  qf.ProjectGridFunction(*GridF);
166 }
167 
169 {
170  if (Q1) { Q1->SetTime(t); }
171  if (Q2) { Q2->SetTime(t); }
172  this->Coefficient::SetTime(t);
173 }
174 
176  const IntegrationPoint &ip)
177 {
178  if (Q2)
179  {
180  return (*Transform2)(Q1->Eval(T, ip, GetTime()),
181  Q2->Eval(T, ip, GetTime()));
182  }
183  else
184  {
185  return (*Transform1)(Q1->Eval(T, ip, GetTime()));
186  }
187 }
188 
190 {
191  if (weight) { weight->SetTime(t); }
192  this->Coefficient::SetTime(t);
193 }
194 
196 {
197  MFEM_VERIFY(vcenter.Size() <= 3,
198  "SetDeltaCenter::Maximum number of dim supported is 3")
199  for (int i = 0; i < vcenter.Size(); i++) { center[i] = vcenter[i]; }
200  sdim = vcenter.Size();
201 }
202 
204 {
205  vcenter.SetSize(sdim);
206  vcenter = center;
207 }
208 
210  const IntegrationPoint &ip)
211 {
212  double w = Scale();
213  return weight ? weight->Eval(T, ip, GetTime())*w : w;
214 }
215 
217 {
218  if (c) { c->SetTime(t); }
219  this->Coefficient::SetTime(t);
220 }
221 
223  const IntegrationRule &ir)
224 {
225  Vector Mi;
226  M.SetSize(vdim, ir.GetNPoints());
227  for (int i = 0; i < ir.GetNPoints(); i++)
228  {
229  M.GetColumnReference(i, Mi);
230  const IntegrationPoint &ip = ir.IntPoint(i);
231  T.SetIntPoint(&ip);
232  Eval(Mi, T, ip);
233  }
234 }
235 
237 {
238  MFEM_VERIFY(vdim == qf.GetVDim(), "Wrong sizes.");
239  QuadratureSpaceBase &qspace = *qf.GetSpace();
240  const int ne = qspace.GetNE();
241  DenseMatrix values;
242  Vector col;
243  for (int iel = 0; iel < ne; ++iel)
244  {
245  qf.GetValues(iel, values);
246  const IntegrationRule &ir = qspace.GetIntRule(iel);
247  ElementTransformation& T = *qspace.GetTransformation(iel);
248  for (int iq = 0; iq < ir.Size(); ++iq)
249  {
250  const IntegrationPoint &ip = ir[iq];
251  T.SetIntPoint(&ip);
252  const int iq_p = qspace.GetPermutedIndex(iel, iq);
253  values.GetColumnReference(iq_p, col);
254  Eval(col, T, ip);
255  }
256  }
257 }
258 
259 void PWVectorCoefficient::InitMap(const Array<int> & attr,
260  const Array<VectorCoefficient*> & coefs)
261 {
262  MFEM_VERIFY(attr.Size() == coefs.Size(),
263  "PWVectorCoefficient: "
264  "Attribute and coefficient arrays have incompatible "
265  "dimensions.");
266 
267  for (int i=0; i<attr.Size(); i++)
268  {
269  if (coefs[i] != NULL)
270  {
271  UpdateCoefficient(attr[i], *coefs[i]);
272  }
273  }
274 }
275 
277 {
278  MFEM_VERIFY(coef.GetVDim() == vdim,
279  "PWVectorCoefficient::UpdateCoefficient: "
280  "VectorCoefficient has incompatible dimension.");
281  pieces[attr] = &coef;
282 }
283 
285 {
287 
288  std::map<int, VectorCoefficient*>::iterator p = pieces.begin();
289  for (; p != pieces.end(); p++)
290  {
291  if (p->second != NULL)
292  {
293  p->second->SetTime(t);
294  }
295  }
296 }
297 
299  const IntegrationPoint &ip)
300 {
301  const int att = T.Attribute;
302  std::map<int, VectorCoefficient*>::const_iterator p = pieces.find(att);
303  if (p != pieces.end())
304  {
305  if ( p->second != NULL)
306  {
307  p->second->Eval(V, T, ip);
308  return;
309  }
310  }
311 
312  V.SetSize(vdim);
313  V = 0.0;
314 }
315 
317  const IntegrationPoint &ip)
318 {
319  double x[3];
320  Vector transip(x, 3);
321 
322  T.Transform(ip, transip);
323 
324  V.SetSize(vdim);
325  if (Function)
326  {
327  Function(transip, V);
328  }
329  else
330  {
331  TDFunction(transip, GetTime(), V);
332  }
333  if (Q)
334  {
335  V *= Q->Eval(T, ip, GetTime());
336  }
337 }
338 
340  : VectorCoefficient(dim), Coeff(dim), ownCoeff(dim)
341 {
342  for (int i = 0; i < dim; i++)
343  {
344  Coeff[i] = NULL;
345  ownCoeff[i] = true;
346  }
347 }
348 
350 {
351  for (int i = 0; i < vdim; i++)
352  {
353  if (Coeff[i]) { Coeff[i]->SetTime(t); }
354  }
356 }
357 
358 void VectorArrayCoefficient::Set(int i, Coefficient *c, bool own)
359 {
360  if (ownCoeff[i]) { delete Coeff[i]; }
361  Coeff[i] = c;
362  ownCoeff[i] = own;
363 }
364 
366 {
367  for (int i = 0; i < vdim; i++)
368  {
369  if (ownCoeff[i]) { delete Coeff[i]; }
370  }
371 }
372 
374  const IntegrationPoint &ip)
375 {
376  V.SetSize(vdim);
377  for (int i = 0; i < vdim; i++)
378  {
379  V(i) = this->Eval(i, T, ip);
380  }
381 }
382 
384  const GridFunction *gf)
385  : VectorCoefficient ((gf) ? gf -> VectorDim() : 0)
386 {
387  GridFunc = gf;
388 }
389 
391 {
392  GridFunc = gf; vdim = (gf) ? gf -> VectorDim() : 0;
393 }
394 
396  const IntegrationPoint &ip)
397 {
398  Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
399  if (T.mesh == gf_mesh)
400  {
401  GridFunc->GetVectorValue(T, ip, V);
402  }
403  else
404  {
405  IntegrationPoint coarse_ip;
406  ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
407  GridFunc->GetVectorValue(*coarse_T, coarse_ip, V);
408  }
409 }
410 
413 {
414  if (T.mesh == GridFunc->FESpace()->GetMesh())
415  {
416  GridFunc->GetVectorValues(T, ir, M);
417  }
418  else
419  {
420  VectorCoefficient::Eval(M, T, ir);
421  }
422 }
423 
425 {
427 }
428 
430  const GridFunction *gf)
431  : VectorCoefficient((gf) ?
432  gf -> FESpace() -> GetMesh() -> SpaceDimension() : 0)
433 {
434  GridFunc = gf;
435 }
436 
438 {
439  GridFunc = gf; vdim = (gf) ?
440  gf -> FESpace() -> GetMesh() -> SpaceDimension() : 0;
441 }
442 
444  const IntegrationPoint &ip)
445 {
446  Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
447  if (T.mesh == gf_mesh)
448  {
449  GridFunc->GetGradient(T, V);
450  }
451  else
452  {
453  IntegrationPoint coarse_ip;
454  ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
455  GridFunc->GetGradient(*coarse_T, V);
456  }
457 }
458 
461 {
462  if (T.mesh == GridFunc->FESpace()->GetMesh())
463  {
464  GridFunc->GetGradients(T, ir, M);
465  }
466  else
467  {
468  VectorCoefficient::Eval(M, T, ir);
469  }
470 }
471 
473  const GridFunction *gf)
474  : VectorCoefficient(0)
475 {
476  SetGridFunction(gf);
477 }
478 
480 {
481  GridFunc = gf; vdim = (gf) ? gf -> CurlDim() : 0;
482 }
483 
485  const IntegrationPoint &ip)
486 {
487  Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
488  if (T.mesh == gf_mesh)
489  {
490  GridFunc->GetCurl(T, V);
491  }
492  else
493  {
494  IntegrationPoint coarse_ip;
495  ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
496  GridFunc->GetCurl(*coarse_T, V);
497  }
498 }
499 
501  const GridFunction *gf) : Coefficient()
502 {
503  GridFunc = gf;
504 }
505 
507  const IntegrationPoint &ip)
508 {
509  Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
510  if (T.mesh == gf_mesh)
511  {
512  return GridFunc->GetDivergence(T);
513  }
514  else
515  {
516  IntegrationPoint coarse_ip;
517  ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
518  return GridFunc->GetDivergence(*coarse_T);
519  }
520 }
521 
523 {
524  d.SetTime(t);
526 }
527 
529 {
530  dir = d_;
531  (*this).vdim = dir.Size();
532 }
533 
536 {
537  V = dir;
538  d.SetTime(GetTime());
539  V *= d.EvalDelta(T, ip);
540 }
541 
543 {
544  if (c) { c->SetTime(t); }
546 }
547 
549  const IntegrationPoint &ip)
550 {
551  V.SetSize(vdim);
552  if (active_attr[T.Attribute-1])
553  {
554  c->SetTime(GetTime());
555  c->Eval(V, T, ip);
556  }
557  else
558  {
559  V = 0.0;
560  }
561 }
562 
565 {
566  if (active_attr[T.Attribute-1])
567  {
568  c->SetTime(GetTime());
569  c->Eval(M, T, ir);
570  }
571  else
572  {
573  M.SetSize(vdim, ir.GetNPoints());
574  M = 0.0;
575  }
576 }
577 
579 {
580  MFEM_VERIFY(qf.GetVDim() == height*width, "Wrong sizes.");
581  QuadratureSpaceBase &qspace = *qf.GetSpace();
582  const int ne = qspace.GetNE();
583  DenseMatrix values, matrix;
584  for (int iel = 0; iel < ne; ++iel)
585  {
586  qf.GetValues(iel, values);
587  const IntegrationRule &ir = qspace.GetIntRule(iel);
588  ElementTransformation& T = *qspace.GetTransformation(iel);
589  for (int iq = 0; iq < ir.Size(); ++iq)
590  {
591  const IntegrationPoint &ip = ir[iq];
592  T.SetIntPoint(&ip);
593  const int iq_p = qspace.GetPermutedIndex(iel, iq);
594  matrix.UseExternalData(&values(0, iq_p), height, width);
595  Eval(matrix, T, ip);
596  if (transpose) { matrix.Transpose(); }
597  }
598  }
599 }
600 
601 void PWMatrixCoefficient::InitMap(const Array<int> & attr,
602  const Array<MatrixCoefficient*> & coefs)
603 {
604  MFEM_VERIFY(attr.Size() == coefs.Size(),
605  "PWMatrixCoefficient: "
606  "Attribute and coefficient arrays have incompatible "
607  "dimensions.");
608 
609  for (int i=0; i<attr.Size(); i++)
610  {
611  if (coefs[i] != NULL)
612  {
613  UpdateCoefficient(attr[i], *coefs[i]);
614  }
615  }
616 }
617 
619 {
620  MFEM_VERIFY(coef.GetHeight() == height,
621  "PWMatrixCoefficient::UpdateCoefficient: "
622  "MatrixCoefficient has incompatible height.");
623  MFEM_VERIFY(coef.GetWidth() == width,
624  "PWMatrixCoefficient::UpdateCoefficient: "
625  "MatrixCoefficient has incompatible width.");
626  if (symmetric)
627  {
628  MFEM_VERIFY(coef.IsSymmetric(),
629  "PWMatrixCoefficient::UpdateCoefficient: "
630  "MatrixCoefficient has incompatible symmetry.");
631  }
632  pieces[attr] = &coef;
633 }
634 
636 {
638 
639  std::map<int, MatrixCoefficient*>::iterator p = pieces.begin();
640  for (; p != pieces.end(); p++)
641  {
642  if (p->second != NULL)
643  {
644  p->second->SetTime(t);
645  }
646  }
647 }
648 
650  const IntegrationPoint &ip)
651 {
652  const int att = T.Attribute;
653  std::map<int, MatrixCoefficient*>::const_iterator p = pieces.find(att);
654  if (p != pieces.end())
655  {
656  if ( p->second != NULL)
657  {
658  p->second->Eval(K, T, ip);
659  return;
660  }
661  }
662 
663  K.SetSize(height, width);
664  K = 0.0;
665 }
666 
668 {
669  if (Q) { Q->SetTime(t); }
671 }
672 
674  const IntegrationPoint &ip)
675 {
676  double x[3];
677  Vector transip(x, 3);
678 
679  T.Transform(ip, transip);
680 
681  K.SetSize(height, width);
682 
683  if (symmetric) // Use SymmFunction (deprecated version)
684  {
685  MFEM_VERIFY(height == width && SymmFunction,
686  "MatrixFunctionCoefficient is not symmetric");
687 
688  Vector Ksym((width * (width + 1)) / 2); // 1x1: 1, 2x2: 3, 3x3: 6
689 
690  SymmFunction(transip, Ksym);
691 
692  // Copy upper triangular values from Ksym to the full matrix K
693  int os = 0;
694  for (int i=0; i<height; ++i)
695  {
696  for (int j=i; j<width; ++j)
697  {
698  const double Kij = Ksym[j - i + os];
699  K(i,j) = Kij;
700  if (j != i) { K(j,i) = Kij; }
701  }
702 
703  os += width - i;
704  }
705  }
706  else
707  {
708  if (Function)
709  {
710  Function(transip, K);
711  }
712  else if (TDFunction)
713  {
714  TDFunction(transip, GetTime(), K);
715  }
716  else
717  {
718  K = mat;
719  }
720  }
721 
722  if (Q)
723  {
724  K *= Q->Eval(T, ip, GetTime());
725  }
726 }
727 
730  const IntegrationPoint &ip)
731 {
732  MFEM_VERIFY(symmetric && height == width && SymmFunction,
733  "MatrixFunctionCoefficient is not symmetric");
734 
735  double x[3];
736  Vector transip(x, 3);
737 
738  T.Transform(ip, transip);
739 
740  K.SetSize((width * (width + 1)) / 2); // 1x1: 1, 2x2: 3, 3x3: 6
741 
742  if (SymmFunction)
743  {
744  SymmFunction(transip, K);
745  }
746 
747  if (Q)
748  {
749  K *= Q->Eval(T, ip, GetTime());
750  }
751 }
752 
754 {
755  const int vdim = qf.GetVDim();
756  MFEM_VERIFY(vdim == height*(height+1)/2, "Wrong sizes.");
757 
758  QuadratureSpaceBase &qspace = *qf.GetSpace();
759  const int ne = qspace.GetNE();
760  DenseMatrix values;
761  DenseSymmetricMatrix matrix;
762  for (int iel = 0; iel < ne; ++iel)
763  {
764  qf.GetValues(iel, values);
765  const IntegrationRule &ir = qspace.GetIntRule(iel);
766  ElementTransformation& T = *qspace.GetTransformation(iel);
767  for (int iq = 0; iq < ir.Size(); ++iq)
768  {
769  const IntegrationPoint &ip = ir[iq];
770  T.SetIntPoint(&ip);
771  matrix.UseExternalData(&values(0, iq), vdim);
772  Eval(matrix, T, ip);
773  }
774  }
775 }
776 
777 
779  const IntegrationPoint &ip)
780 {
781  mat.SetSize(height);
782  Eval(mat, T, ip);
783  for (int j = 0; j < width; ++j)
784  {
785  for (int i = 0; i < height; ++ i)
786  {
787  K(i, j) = mat(i, j);
788  }
789  }
790 }
791 
793 {
794  if (Q) { Q->SetTime(t); }
796 }
797 
800  const IntegrationPoint &ip)
801 {
802  double x[3];
803  Vector transip(x, 3);
804 
805  T.Transform(ip, transip);
806 
807  K.SetSize(height);
808 
809  if (Function)
810  {
811  Function(transip, K);
812  }
813  else if (TDFunction)
814  {
815  TDFunction(transip, GetTime(), K);
816  }
817  else
818  {
819  K = mat;
820  }
821 
822  if (Q)
823  {
824  K *= Q->Eval(T, ip, GetTime());
825  }
826 }
827 
830 {
831  Coeff.SetSize(height*width);
832  ownCoeff.SetSize(height*width);
833  for (int i = 0; i < (height*width); i++)
834  {
835  Coeff[i] = NULL;
836  ownCoeff[i] = true;
837  }
838 }
839 
841 {
842  for (int i=0; i < height*width; i++)
843  {
844  if (Coeff[i]) { Coeff[i]->SetTime(t); }
845  }
847 }
848 
849 void MatrixArrayCoefficient::Set(int i, int j, Coefficient * c, bool own)
850 {
851  if (ownCoeff[i*width+j]) { delete Coeff[i*width+j]; }
852  Coeff[i*width+j] = c;
853  ownCoeff[i*width+j] = own;
854 }
855 
857 {
858  for (int i=0; i < height*width; i++)
859  {
860  if (ownCoeff[i]) { delete Coeff[i]; }
861  }
862 }
863 
865  const IntegrationPoint &ip)
866 {
867  K.SetSize(height, width);
868  for (int i = 0; i < height; i++)
869  {
870  for (int j = 0; j < width; j++)
871  {
872  K(i,j) = this->Eval(i, j, T, ip);
873  }
874  }
875 }
876 
878 {
879  if (c) { c->SetTime(t); }
881 }
882 
884  const IntegrationPoint &ip)
885 {
886  if (active_attr[T.Attribute-1])
887  {
888  c->SetTime(GetTime());
889  c->Eval(K, T, ip);
890  }
891  else
892  {
893  K.SetSize(height, width);
894  K = 0.0;
895  }
896 }
897 
899 {
900  if (a) { a->SetTime(t); }
901  if (b) { b->SetTime(t); }
902  this->Coefficient::SetTime(t);
903 }
904 
906 {
907  if (a) { a->SetTime(t); }
908  if (b) { b->SetTime(t); }
909  this->Coefficient::SetTime(t);
910 }
911 
913 {
914  if (a) { a->SetTime(t); }
915  if (b) { b->SetTime(t); }
916  this->Coefficient::SetTime(t);
917 }
918 
920 {
921  if (a) { a->SetTime(t); }
922  this->Coefficient::SetTime(t);
923 }
924 
927  : a(&A), b(&B), va(A.GetVDim()), vb(B.GetVDim())
928 {
929  MFEM_ASSERT(A.GetVDim() == B.GetVDim(),
930  "InnerProductCoefficient: "
931  "Arguments have incompatible dimensions.");
932 }
933 
935 {
936  if (a) { a->SetTime(t); }
937  if (b) { b->SetTime(t); }
938  this->Coefficient::SetTime(t);
939 }
940 
942  const IntegrationPoint &ip)
943 {
944  a->Eval(va, T, ip);
945  b->Eval(vb, T, ip);
946  return va * vb;
947 }
948 
951  : a(&A), b(&B), va(A.GetVDim()), vb(B.GetVDim())
952 {
953  MFEM_ASSERT(A.GetVDim() == 2 && B.GetVDim() == 2,
954  "VectorRotProductCoefficient: "
955  "Arguments must have dimension equal to two.");
956 }
957 
959 {
960  if (a) { a->SetTime(t); }
961  if (b) { b->SetTime(t); }
962  this->Coefficient::SetTime(t);
963 }
964 
966  const IntegrationPoint &ip)
967 {
968  a->Eval(va, T, ip);
969  b->Eval(vb, T, ip);
970  return va[0] * vb[1] - va[1] * vb[0];
971 }
972 
974  : a(&A), ma(A.GetHeight(), A.GetWidth())
975 {
976  MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
977  "DeterminantCoefficient: "
978  "Argument must be a square matrix.");
979 }
980 
982 {
983  if (a) { a->SetTime(t); }
984  this->Coefficient::SetTime(t);
985 }
986 
988  const IntegrationPoint &ip)
989 {
990  a->Eval(ma, T, ip);
991  return ma.Det();
992 }
993 
996  ACoef(NULL), BCoef(NULL),
997  A(dim), B(dim),
998  alphaCoef(NULL), betaCoef(NULL),
999  alpha(1.0), beta(1.0)
1000 {
1001  A = 0.0; B = 0.0;
1002 }
1003 
1005  VectorCoefficient &B_,
1006  double alpha_, double beta_)
1007  : VectorCoefficient(A_.GetVDim()),
1008  ACoef(&A_), BCoef(&B_),
1009  A(A_.GetVDim()), B(A_.GetVDim()),
1010  alphaCoef(NULL), betaCoef(NULL),
1011  alpha(alpha_), beta(beta_)
1012 {
1013  MFEM_ASSERT(A_.GetVDim() == B_.GetVDim(),
1014  "VectorSumCoefficient: "
1015  "Arguments must have the same dimension.");
1016 }
1017 
1019  VectorCoefficient &B_,
1020  Coefficient &alpha_,
1021  Coefficient &beta_)
1022  : VectorCoefficient(A_.GetVDim()),
1023  ACoef(&A_), BCoef(&B_),
1024  A(A_.GetVDim()),
1025  B(A_.GetVDim()),
1026  alphaCoef(&alpha_),
1027  betaCoef(&beta_),
1028  alpha(0.0), beta(0.0)
1029 {
1030  MFEM_ASSERT(A_.GetVDim() == B_.GetVDim(),
1031  "VectorSumCoefficient: "
1032  "Arguments must have the same dimension.");
1033 }
1034 
1036 {
1037  if (ACoef) { ACoef->SetTime(t); }
1038  if (BCoef) { BCoef->SetTime(t); }
1039  if (alphaCoef) { alphaCoef->SetTime(t); }
1040  if (betaCoef) { betaCoef->SetTime(t); }
1041  this->VectorCoefficient::SetTime(t);
1042 }
1043 
1045  const IntegrationPoint &ip)
1046 {
1047  V.SetSize(A.Size());
1048  if ( ACoef) { ACoef->Eval(A, T, ip); }
1049  if ( BCoef) { BCoef->Eval(B, T, ip); }
1050  if (alphaCoef) { alpha = alphaCoef->Eval(T, ip); }
1051  if ( betaCoef) { beta = betaCoef->Eval(T, ip); }
1052  add(alpha, A, beta, B, V);
1053 }
1054 
1056  double A,
1057  VectorCoefficient &B)
1058  : VectorCoefficient(B.GetVDim()), aConst(A), a(NULL), b(&B)
1059 {}
1060 
1062  Coefficient &A,
1063  VectorCoefficient &B)
1064  : VectorCoefficient(B.GetVDim()), aConst(0.0), a(&A), b(&B)
1065 {}
1066 
1068 {
1069  if (a) { a->SetTime(t); }
1070  if (b) { b->SetTime(t); }
1071  this->VectorCoefficient::SetTime(t);
1072 }
1073 
1075  const IntegrationPoint &ip)
1076 {
1077  double sa = (a == NULL) ? aConst : a->Eval(T, ip);
1078  b->Eval(V, T, ip);
1079  V *= sa;
1080 }
1081 
1083  double tol_)
1084  : VectorCoefficient(A.GetVDim()), a(&A), tol(tol_)
1085 {}
1086 
1088 {
1089  if (a) { a->SetTime(t); }
1090  this->VectorCoefficient::SetTime(t);
1091 }
1092 
1094  const IntegrationPoint &ip)
1095 {
1096  a->Eval(V, T, ip);
1097  double nv = V.Norml2();
1098  V *= (nv > tol) ? (1.0/nv) : 0.0;
1099 }
1100 
1102  VectorCoefficient &A,
1103  VectorCoefficient &B)
1104  : VectorCoefficient(3), a(&A), b(&B), va(A.GetVDim()), vb(B.GetVDim())
1105 {
1106  MFEM_ASSERT(A.GetVDim() == 3 && B.GetVDim() == 3,
1107  "VectorCrossProductCoefficient: "
1108  "Arguments must have dimension equal to three.");
1109 }
1110 
1112 {
1113  if (a) { a->SetTime(t); }
1114  if (b) { b->SetTime(t); }
1115  this->VectorCoefficient::SetTime(t);
1116 }
1117 
1119  const IntegrationPoint &ip)
1120 {
1121  a->Eval(va, T, ip);
1122  b->Eval(vb, T, ip);
1123  V.SetSize(3);
1124  V[0] = va[1] * vb[2] - va[2] * vb[1];
1125  V[1] = va[2] * vb[0] - va[0] * vb[2];
1126  V[2] = va[0] * vb[1] - va[1] * vb[0];
1127 }
1128 
1131  : VectorCoefficient(A.GetHeight()), a(&A), b(&B),
1132  ma(A.GetHeight(), A.GetWidth()), vb(B.GetVDim())
1133 {
1134  MFEM_ASSERT(A.GetWidth() == B.GetVDim(),
1135  "MatrixVectorProductCoefficient: "
1136  "Arguments have incompatible dimensions.");
1137 }
1138 
1140 {
1141  if (a) { a->SetTime(t); }
1142  if (b) { b->SetTime(t); }
1143  this->VectorCoefficient::SetTime(t);
1144 }
1145 
1147  const IntegrationPoint &ip)
1148 {
1149  a->Eval(ma, T, ip);
1150  b->Eval(vb, T, ip);
1151  V.SetSize(vdim);
1152  ma.Mult(vb, V);
1153 }
1154 
1156  const IntegrationPoint &ip)
1157 {
1158  M.SetSize(dim);
1159  M = 0.0;
1160  for (int d=0; d<dim; d++) { M(d,d) = 1.0; }
1161 }
1162 
1164  MatrixCoefficient &B,
1165  double alpha_, double beta_)
1166  : MatrixCoefficient(A.GetHeight(), A.GetWidth()),
1167  a(&A), b(&B), alpha(alpha_), beta(beta_),
1168  ma(A.GetHeight(), A.GetWidth())
1169 {
1170  MFEM_ASSERT(A.GetHeight() == B.GetHeight() && A.GetWidth() == B.GetWidth(),
1171  "MatrixSumCoefficient: "
1172  "Arguments must have the same dimensions.");
1173 }
1174 
1176 {
1177  if (a) { a->SetTime(t); }
1178  if (b) { b->SetTime(t); }
1179  this->MatrixCoefficient::SetTime(t);
1180 }
1181 
1183  const IntegrationPoint &ip)
1184 {
1185  b->Eval(M, T, ip);
1186  if ( beta != 1.0 ) { M *= beta; }
1187  a->Eval(ma, T, ip);
1188  M.Add(alpha, ma);
1189 }
1190 
1192  MatrixCoefficient &B)
1193  : MatrixCoefficient(A.GetHeight(), B.GetWidth()),
1194  a(&A), b(&B),
1195  ma(A.GetHeight(), A.GetWidth()),
1196  mb(B.GetHeight(), B.GetWidth())
1197 {
1198  MFEM_ASSERT(A.GetWidth() == B.GetHeight(),
1199  "MatrixProductCoefficient: "
1200  "Arguments must have compatible dimensions.");
1201 }
1202 
1204  const IntegrationPoint &ip)
1205 {
1206  a->Eval(ma, T, ip);
1207  b->Eval(mb, T, ip);
1208  Mult(ma, mb, M);
1209 }
1210 
1212  double A,
1213  MatrixCoefficient &B)
1214  : MatrixCoefficient(B.GetHeight(), B.GetWidth()), aConst(A), a(NULL), b(&B)
1215 {}
1216 
1218  Coefficient &A,
1219  MatrixCoefficient &B)
1220  : MatrixCoefficient(B.GetHeight(), B.GetWidth()), aConst(0.0), a(&A), b(&B)
1221 {}
1222 
1224 {
1225  if (a) { a->SetTime(t); }
1226  if (b) { b->SetTime(t); }
1227  this->MatrixCoefficient::SetTime(t);
1228 }
1229 
1232  const IntegrationPoint &ip)
1233 {
1234  double sa = (a == NULL) ? aConst : a->Eval(T, ip);
1235  b->Eval(M, T, ip);
1236  M *= sa;
1237 }
1238 
1240  : MatrixCoefficient(A.GetWidth(), A.GetHeight()), a(&A)
1241 {}
1242 
1244 {
1245  if (a) { a->SetTime(t); }
1246  this->MatrixCoefficient::SetTime(t);
1247 }
1248 
1251  const IntegrationPoint &ip)
1252 {
1253  a->Eval(M, T, ip);
1254  M.Transpose();
1255 }
1256 
1258  : MatrixCoefficient(A.GetHeight(), A.GetWidth()), a(&A)
1259 {
1260  MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
1261  "InverseMatrixCoefficient: "
1262  "Argument must be a square matrix.");
1263 }
1264 
1266 {
1267  if (a) { a->SetTime(t); }
1268  this->MatrixCoefficient::SetTime(t);
1269 }
1270 
1273  const IntegrationPoint &ip)
1274 {
1275  a->Eval(M, T, ip);
1276  M.Invert();
1277 }
1278 
1280  VectorCoefficient &B)
1281  : MatrixCoefficient(A.GetVDim(), B.GetVDim()), a(&A), b(&B),
1282  va(A.GetVDim()), vb(B.GetVDim())
1283 {}
1284 
1286 {
1287  if (a) { a->SetTime(t); }
1288  if (b) { b->SetTime(t); }
1289  this->MatrixCoefficient::SetTime(t);
1290 }
1291 
1293  const IntegrationPoint &ip)
1294 {
1295  a->Eval(va, T, ip);
1296  b->Eval(vb, T, ip);
1297  M.SetSize(va.Size(), vb.Size());
1298  for (int i=0; i<va.Size(); i++)
1299  {
1300  for (int j=0; j<vb.Size(); j++)
1301  {
1302  M(i, j) = va[i] * vb[j];
1303  }
1304  }
1305 }
1306 
1308  : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(A), a(NULL), k(&K),
1309  vk(K.GetVDim())
1310 {}
1311 
1313  VectorCoefficient &K)
1314  : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(0.0), a(&A), k(&K),
1315  vk(K.GetVDim())
1316 {}
1317 
1319 {
1320  if (a) { a->SetTime(t); }
1321  if (k) { k->SetTime(t); }
1322  this->MatrixCoefficient::SetTime(t);
1323 }
1324 
1326  const IntegrationPoint &ip)
1327 {
1328  k->Eval(vk, T, ip);
1329  M.SetSize(vk.Size(), vk.Size());
1330  M = 0.0;
1331  double k2 = vk*vk;
1332  for (int i=0; i<vk.Size(); i++)
1333  {
1334  M(i, i) = k2;
1335  for (int j=0; j<vk.Size(); j++)
1336  {
1337  M(i, j) -= vk[i] * vk[j];
1338  }
1339  }
1340  M *= ((a == NULL ) ? aConst : a->Eval(T, ip) );
1341 }
1342 
1343 double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh,
1344  const IntegrationRule *irs[])
1345 {
1346  double norm = 0.0;
1348 
1349  for (int i = 0; i < mesh.GetNE(); i++)
1350  {
1351  tr = mesh.GetElementTransformation(i);
1352  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
1353  for (int j = 0; j < ir.GetNPoints(); j++)
1354  {
1355  const IntegrationPoint &ip = ir.IntPoint(j);
1356  tr->SetIntPoint(&ip);
1357  double val = fabs(coeff.Eval(*tr, ip));
1358  if (p < infinity())
1359  {
1360  norm += ip.weight * tr->Weight() * pow(val, p);
1361  }
1362  else
1363  {
1364  if (norm < val)
1365  {
1366  norm = val;
1367  }
1368  }
1369  }
1370  }
1371  return norm;
1372 }
1373 
1374 double LpNormLoop(double p, VectorCoefficient &coeff, Mesh &mesh,
1375  const IntegrationRule *irs[])
1376 {
1377  double norm = 0.0;
1379  int vdim = coeff.GetVDim();
1380  Vector vval(vdim);
1381  double val;
1382 
1383  for (int i = 0; i < mesh.GetNE(); i++)
1384  {
1385  tr = mesh.GetElementTransformation(i);
1386  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
1387  for (int j = 0; j < ir.GetNPoints(); j++)
1388  {
1389  const IntegrationPoint &ip = ir.IntPoint(j);
1390  tr->SetIntPoint(&ip);
1391  coeff.Eval(vval, *tr, ip);
1392  if (p < infinity())
1393  {
1394  for (int idim(0); idim < vdim; ++idim)
1395  {
1396  norm += ip.weight * tr->Weight() * pow(fabs( vval(idim) ), p);
1397  }
1398  }
1399  else
1400  {
1401  for (int idim(0); idim < vdim; ++idim)
1402  {
1403  val = fabs(vval(idim));
1404  if (norm < val)
1405  {
1406  norm = val;
1407  }
1408  }
1409  }
1410  }
1411  }
1412 
1413  return norm;
1414 }
1415 
1416 double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh,
1417  const IntegrationRule *irs[])
1418 {
1419  double norm = LpNormLoop(p, coeff, mesh, irs);
1420 
1421  if (p < infinity())
1422  {
1423  // negative quadrature weights may cause norm to be negative
1424  if (norm < 0.0)
1425  {
1426  norm = -pow(-norm, 1.0/p);
1427  }
1428  else
1429  {
1430  norm = pow(norm, 1.0/p);
1431  }
1432  }
1433 
1434  return norm;
1435 }
1436 
1437 double ComputeLpNorm(double p, VectorCoefficient &coeff, Mesh &mesh,
1438  const IntegrationRule *irs[])
1439 {
1440  double norm = LpNormLoop(p, coeff, mesh, irs);
1441 
1442  if (p < infinity())
1443  {
1444  // negative quadrature weights may cause norm to be negative
1445  if (norm < 0.0)
1446  {
1447  norm = -pow(-norm, 1.0/p);
1448  }
1449  else
1450  {
1451  norm = pow(norm, 1.0/p);
1452  }
1453  }
1454 
1455  return norm;
1456 }
1457 
1458 #ifdef MFEM_USE_MPI
1459 double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh,
1460  const IntegrationRule *irs[])
1461 {
1462  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
1463  double glob_norm = 0;
1464 
1465  MPI_Comm comm = pmesh.GetComm();
1466 
1467  if (p < infinity())
1468  {
1469  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
1470 
1471  // negative quadrature weights may cause norm to be negative
1472  if (glob_norm < 0.0)
1473  {
1474  glob_norm = -pow(-glob_norm, 1.0/p);
1475  }
1476  else
1477  {
1478  glob_norm = pow(glob_norm, 1.0/p);
1479  }
1480  }
1481  else
1482  {
1483  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
1484  }
1485 
1486  return glob_norm;
1487 }
1488 
1489 double ComputeGlobalLpNorm(double p, VectorCoefficient &coeff, ParMesh &pmesh,
1490  const IntegrationRule *irs[])
1491 {
1492  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
1493  double glob_norm = 0;
1494 
1495  MPI_Comm comm = pmesh.GetComm();
1496 
1497  if (p < infinity())
1498  {
1499  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
1500 
1501  // negative quadrature weights may cause norm to be negative
1502  if (glob_norm < 0.0)
1503  {
1504  glob_norm = -pow(-glob_norm, 1.0/p);
1505  }
1506  else
1507  {
1508  glob_norm = pow(glob_norm, 1.0/p);
1509  }
1510  }
1511  else
1512  {
1513  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
1514  }
1515 
1516  return glob_norm;
1517 }
1518 #endif
1519 
1521  QuadratureFunction &qf)
1522  : VectorCoefficient(qf.GetVDim()), QuadF(qf), index(0) { }
1523 
1525 {
1526  MFEM_VERIFY(index_ >= 0, "Index must be >= 0");
1527  MFEM_VERIFY(index_ < QuadF.GetVDim(),
1528  "Index must be < QuadratureFunction length");
1529  index = index_;
1530 
1531  MFEM_VERIFY(length_ > 0, "Length must be > 0");
1532  MFEM_VERIFY(length_ <= QuadF.GetVDim() - index,
1533  "Length must be <= (QuadratureFunction length - index)");
1534 
1535  vdim = length_;
1536 }
1537 
1540  const IntegrationPoint &ip)
1541 {
1542  QuadF.HostRead();
1543 
1544  if (index == 0 && vdim == QuadF.GetVDim())
1545  {
1546  QuadF.GetValues(T.ElementNo, ip.index, V);
1547  }
1548  else
1549  {
1550  Vector temp;
1551  QuadF.GetValues(T.ElementNo, ip.index, temp);
1552  V.SetSize(vdim);
1553  for (int i = 0; i < vdim; i++)
1554  {
1555  V(i) = temp(index + i);
1556  }
1557  }
1558 
1559  return;
1560 }
1561 
1563 {
1564  qf = QuadF;
1565 }
1566 
1568  QuadratureFunction &qf) : QuadF(qf)
1569 {
1570  MFEM_VERIFY(qf.GetVDim() == 1, "QuadratureFunction's vdim must be 1");
1571 }
1572 
1574  const IntegrationPoint &ip)
1575 {
1576  QuadF.HostRead();
1577  Vector temp(1);
1578  QuadF.GetValues(T.ElementNo, ip.index, temp);
1579  return temp[0];
1580 }
1581 
1583 {
1584  qf = QuadF;
1585 }
1586 
1587 
1589  QuadratureSpaceBase &qs_, CoefficientStorage storage_)
1590  : Vector(), storage(storage_), vdim(0), qs(qs_), qf(NULL)
1591 {
1592  UseDevice(true);
1593 }
1594 
1596  QuadratureSpaceBase &qs_,
1597  CoefficientStorage storage_)
1598  : CoefficientVector(qs_, storage_)
1599 {
1600  if (coeff == NULL)
1601  {
1602  SetConstant(1.0);
1603  }
1604  else
1605  {
1606  Project(*coeff);
1607  }
1608 }
1609 
1611  QuadratureSpaceBase &qs_,
1612  CoefficientStorage storage_)
1613  : CoefficientVector(qs_, storage_)
1614 {
1615  Project(coeff);
1616 }
1617 
1619  QuadratureSpaceBase &qs_,
1620  CoefficientStorage storage_)
1621  : CoefficientVector(qs_, storage_)
1622 {
1623  Project(coeff);
1624 }
1625 
1627  QuadratureSpaceBase &qs_,
1628  CoefficientStorage storage_)
1629  : CoefficientVector(qs_, storage_)
1630 {
1631  Project(coeff);
1632 }
1633 
1635 {
1636  vdim = 1;
1637  if (auto *const_coeff = dynamic_cast<ConstantCoefficient*>(&coeff))
1638  {
1639  SetConstant(const_coeff->constant);
1640  }
1641  else if (auto *qf_coeff = dynamic_cast<QuadratureFunctionCoefficient*>(&coeff))
1642  {
1643  MakeRef(qf_coeff->GetQuadFunction());
1644  }
1645  else
1646  {
1647  if (qf == nullptr) { qf = new QuadratureFunction(qs); }
1648  qf->SetVDim(1);
1649  coeff.Project(*qf);
1650  Vector::MakeRef(*qf, 0, qf->Size());
1651  }
1652 }
1653 
1655 {
1656  vdim = coeff.GetVDim();
1657  if (auto *const_coeff = dynamic_cast<VectorConstantCoefficient*>(&coeff))
1658  {
1659  SetConstant(const_coeff->GetVec());
1660  }
1661  else if (auto *qf_coeff =
1662  dynamic_cast<VectorQuadratureFunctionCoefficient*>(&coeff))
1663  {
1664  MakeRef(qf_coeff->GetQuadFunction());
1665  }
1666  else
1667  {
1668  if (qf == nullptr) { qf = new QuadratureFunction(qs, vdim); }
1669  qf->SetVDim(vdim);
1670  coeff.Project(*qf);
1671  Vector::MakeRef(*qf, 0, qf->Size());
1672  }
1673 }
1674 
1675 void CoefficientVector::Project(MatrixCoefficient &coeff, bool transpose)
1676 {
1677  if (auto *const_coeff = dynamic_cast<MatrixConstantCoefficient*>(&coeff))
1678  {
1679  SetConstant(const_coeff->GetMatrix());
1680  }
1681  else if (auto *const_sym_coeff =
1682  dynamic_cast<SymmetricMatrixConstantCoefficient*>(&coeff))
1683  {
1684  SetConstant(const_sym_coeff->GetMatrix());
1685  }
1686  else
1687  {
1688  auto *sym_coeff = dynamic_cast<SymmetricMatrixCoefficient*>(&coeff);
1689  const bool sym = sym_coeff && (storage & CoefficientStorage::SYMMETRIC);
1690  const int height = coeff.GetHeight();
1691  const int width = coeff.GetWidth();
1692  vdim = sym ? height*(height + 1)/2 : width*height;
1693 
1694  if (qf == nullptr) { qf = new QuadratureFunction(qs, vdim); }
1695  qf->SetVDim(vdim);
1696  if (sym) { sym_coeff->ProjectSymmetric(*qf); }
1697  else { coeff.Project(*qf, transpose); }
1698  Vector::MakeRef(*qf, 0, qf->Size());
1699  }
1700 }
1701 
1703 {
1704  Project(coeff, true);
1705 }
1706 
1708 {
1709  vdim = qf_.GetVDim();
1710  const QuadratureSpaceBase *qs2 = qf_.GetSpace();
1711  MFEM_CONTRACT_VAR(qs2); // qs2 used only for asserts
1712  MFEM_VERIFY(qs2 != NULL, "Invalid QuadratureSpace.")
1713  MFEM_VERIFY(qs2->GetMesh() == qs.GetMesh(), "Meshes differ.");
1714  MFEM_VERIFY(qs2->GetOrder() == qs.GetOrder(), "Orders differ.");
1715  Vector::MakeRef(const_cast<QuadratureFunction&>(qf_), 0, qf_.Size());
1716 }
1717 
1718 void CoefficientVector::SetConstant(double constant)
1719 {
1720  const int nq = (storage & CoefficientStorage::CONSTANTS) ? 1 : qs.GetSize();
1721  vdim = 1;
1722  SetSize(nq);
1723  Vector::operator=(constant);
1724 }
1725 
1727 {
1728  const int nq = (storage & CoefficientStorage::CONSTANTS) ? 1 : qs.GetSize();
1729  vdim = constant.Size();
1730  SetSize(nq*vdim);
1731  for (int iq = 0; iq < nq; ++iq)
1732  {
1733  for (int vd = 0; vd<vdim; ++vd)
1734  {
1735  (*this)[vd + iq*vdim] = constant[vd];
1736  }
1737  }
1738 }
1739 
1741 {
1742  const int nq = (storage & CoefficientStorage::CONSTANTS) ? 1 : qs.GetSize();
1743  const int width = constant.Width();
1744  const int height = constant.Height();
1745  vdim = width*height;
1746  SetSize(nq*vdim);
1747  for (int iq = 0; iq < nq; ++iq)
1748  {
1749  for (int j = 0; j < width; ++j)
1750  {
1751  for (int i = 0; i < height; ++i)
1752  {
1753  (*this)[i + j*height + iq*vdim] = constant(i, j);
1754  }
1755  }
1756  }
1757 }
1758 
1760 {
1761  const int nq = (storage & CoefficientStorage::CONSTANTS) ? 1 : qs.GetSize();
1762  const int height = constant.Height();
1763  const bool sym = storage & CoefficientStorage::SYMMETRIC;
1764  vdim = sym ? height*(height + 1)/2 : height*height;
1765  SetSize(nq*vdim);
1766  for (int iq = 0; iq < nq; ++iq)
1767  {
1768  for (int vd = 0; vd < vdim; ++vd)
1769  {
1770  const double value = sym ? constant.GetData()[vd] : constant(vd % height,
1771  vd / height);
1772  (*this)[vd + iq*vdim] = value;
1773  }
1774  }
1775 }
1776 
1777 int CoefficientVector::GetVDim() const { return vdim; }
1778 
1780 {
1781  delete qf;
1782 }
1783 
1784 }
int GetHeight() const
Get the height of the matrix.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient at ip.
CoefficientStorage
Flags that determine what storage optimizations to use in CoefficientVector.
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
int GetNE() const
Return the number of entities.
Definition: qspace.hpp:60
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient at ip.
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void ProjectGridFunction(const GridFunction &gf)
Evaluate a grid function at each quadrature point.
Definition: qfunction.cpp:56
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:479
void Set(int i, int j, Coefficient *c, bool own=true)
Set the coefficient located at (i,j) in the matrix. By default by default this will take ownership of...
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Base class for vector Coefficients that optionally depend on time and space.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector curl coefficient at ip.
TransposeMatrixCoefficient(MatrixCoefficient &A)
Construct with the matrix coefficient. Result is .
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 ...
VectorCrossProductCoefficient(VectorCoefficient &A, VectorCoefficient &B)
Construct with the two coefficients. Result is A x B.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
NormalizedVectorCoefficient(VectorCoefficient &A, double tol=1e-6)
Return a vector normalized to a length of one.
int vdim
Number of values per quadrature point.
void SetSize(int s)
Change the size of the DenseSymmetricMatrix to s x s.
Definition: symmat.cpp:32
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
MatrixVectorProductCoefficient(MatrixCoefficient &A, VectorCoefficient &B)
Constructor with two coefficients. Result is A*B.
double Eval(int i, ElementTransformation &T, const IntegrationPoint &ip)
void SetDeltaCenter(const Vector &center)
Set the center location of the delta function.
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:452
void SetTime(double t)
Set the time for internally stored coefficients.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
void SetTime(double t)
Set the time for internally stored coefficients.
virtual int GetPermutedIndex(int idx, int iq) const =0
Returns the permuted index of the iq quadrature point in entity idx.
void SetTime(double t)
Set the time for internally stored coefficients.
double GetTime()
Get the time for time dependent coefficients.
void GetGradient(ElementTransformation &tr, Vector &grad) const
Gradient of a scalar function at a quadrature point.
Definition: gridfunc.cpp:1646
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
void SetTime(double t)
Set the time for internally stored coefficients.
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
void SetTime(double t)
Set the time for internally stored coefficients.
double Det() const
Definition: densemat.cpp:488
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
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient at ip.
CurlGridFunctionCoefficient(const GridFunction *gf)
Construct the coefficient with a vector grid function gf. The grid function is not owned by the coeff...
void UpdateCoefficient(int attr, VectorCoefficient &coef)
Replace a single Coefficient for a particular attribute.
GradientGridFunctionCoefficient(const GridFunction *gf)
Construct the coefficient with a scalar grid function gf. The grid function is not owned by the coeff...
CrossCrossCoefficient(double A, VectorCoefficient &K)
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
STL namespace.
void SetConstant(double constant)
Set this vector to the given constant.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
virtual void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
virtual void SetTime(double t)
Set the time for time dependent coefficients.
Definition: coefficient.hpp:50
void SetTime(double t)
Set the time for internally stored coefficients.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
int GetWidth() const
Get the width of the matrix.
Class to represent a coefficient evaluated at quadrature points.
virtual ElementTransformation * GetTransformation(int idx)=0
Get the (element or face) transformation of entity idx.
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:314
DeterminantCoefficient(MatrixCoefficient &A)
Construct with the matrix.
ScalarMatrixProductCoefficient(double A, MatrixCoefficient &B)
Constructor with one coefficient. Result is A*B.
virtual void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
void Project(Coefficient &coeff)
Evaluate the given Coefficient at the quadrature points defined by qs.
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:124
virtual void SetTime(double t)
Set the time for time dependent coefficients.
double Eval(int i, int j, ElementTransformation &T, const IntegrationPoint &ip)
InnerProductCoefficient(VectorCoefficient &A, VectorCoefficient &B)
Construct with the two vector coefficients. Result is .
void SetTime(double t)
Set the time for internally stored coefficients.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void SetPointMat(const DenseMatrix &pm)
Set the underlying point matrix describing the transformation.
Definition: eltrans.hpp:401
void SetTime(double t)
Set the time for internally stored coefficients.
void SetTime(double t)
Set the time for internally stored coefficients.
void SetTime(double t)
Set the time for internally stored coefficients.
void SetTime(double t)
Set the time for internally stored coefficients.
void GetCurl(ElementTransformation &tr, Vector &curl) const
Definition: gridfunc.cpp:1550
double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
Compute the Lp norm of a function f. .
void SetTime(double t)
Set the time for internally stored coefficients.
void SetTime(double t)
Set the time for internally stored coefficients.
Store constants using only vdim entries.
DenseTensor point_matrices[Geometry::NumGeom]
Definition: ncmesh.hpp:77
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
double GetTime()
Get the time for time dependent coefficients.
ElementTransformation * RefinedToCoarse(Mesh &coarse_mesh, const ElementTransformation &T, const IntegrationPoint &ip, IntegrationPoint &coarse_ip)
Definition: coefficient.cpp:27
double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
QuadratureFunctionCoefficient(QuadratureFunction &qf)
Constructor with a quadrature function as input.
void Add(const double c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Definition: densemat.cpp:581
double b
Definition: lissajous.cpp:42
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient.
int GetVDim() const
Get the vector dimension.
Definition: qfunction.hpp:74
void UseExternalData(double *d, int s)
Change the data array and the size of the DenseSymmetricMatrix.
Definition: symmat.hpp:48
OuterProductCoefficient(VectorCoefficient &A, VectorCoefficient &B)
Construct with two vector coefficients. Result is .
void SetGridFunction(const GridFunction *gf)
Set the vector grid function.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:679
void Set(int i, Coefficient *c, bool own=true)
Sets coefficient in the vector.
Array< Embedding > embeddings
Fine element positions in their parents.
Definition: ncmesh.hpp:73
Mesh * GetMesh() const
Returns the mesh.
Definition: qspace.hpp:63
void MakeRef(const QuadratureFunction &qf_)
Make this vector a reference to the given QuadratureFunction.
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
void UpdateCoefficient(int attr, MatrixCoefficient &coef)
Replace a single coefficient for a particular attribute.
void SetVDim(int vdim_)
Set the vector dimension, updating the size by calling Vector::SetSize().
Definition: qfunction.hpp:77
virtual ~VectorArrayCoefficient()
Destroys vector coefficient.
void SetTime(double t)
Set the time for internally stored coefficients.
void SetTime(double t)
Set the time for internally stored coefficients.
double * GetData() const
Returns the matrix data array.
Definition: symmat.hpp:80
void SetTime(double t)
Set the time for internally stored coefficients.
Store the triangular part of symmetric matrices.
double GetDivergence(ElementTransformation &tr) const
Definition: gridfunc.cpp:1460
int GetVDim()
Returns dimension of the vector.
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the determinant coefficient at ip.
void SetTime(double t)
Set the time for internally stored coefficients.
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
void SetComponent(int index_, int length_)
Abstract base class for QuadratureSpace and FaceQuadratureSpace.
Definition: qspace.hpp:24
void GetDeltaCenter(Vector &center)
Write the center of the delta function into center.
int GetVDim() const
Return the number of values per quadrature point.
virtual void ProjectSymmetric(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
void SetGridFunction(const GridFunction *gf)
Set the scalar grid function.
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1420
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition: eltrans.cpp:383
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
VectorQuadratureFunctionCoefficient(QuadratureFunction &qf)
Constructor with a quadrature function as input.
virtual void SetTime(double t)
Set the time for time dependent coefficients.
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
virtual void SetTime(double t)
Set the time for time dependent coefficients.
Base class for Matrix Coefficients that optionally depend on time and space.
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:546
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
CoefficientStorage storage
Storage optimizations (see CoefficientStorage).
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient.
Definition: coefficient.cpp:76
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
virtual void Eval(DenseSymmetricMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result as ...
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient at ip.
const IntegrationRule & GetIntRule(int idx) const
Return the IntegrationRule associated with entity idx.
Definition: qspace.hpp:72
virtual void EvalDelta(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Return the specified direction vector multiplied by the value returned by DeltaCoefficient::EvalDelta...
MatrixProductCoefficient(MatrixCoefficient &A, MatrixCoefficient &B)
Construct with the two coefficients. Result is A * B.
double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
Compute the global Lp norm of a function f. .
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void GetValues(int idx, Vector &values)
Return all values associated with mesh element idx in a Vector.
Definition: qfunction.hpp:205
virtual double EvalDelta(ElementTransformation &T, const IntegrationPoint &ip)
The value of the function assuming we are evaluating at the delta center.
const CoarseFineTransformations & GetRefinementTransforms()
Definition: mesh.cpp:9990
virtual void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
virtual void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
double a
Definition: lissajous.cpp:41
MatrixArrayCoefficient(int dim)
Construct a coefficient matrix of dimensions dim * dim. The actual coefficients still need to be adde...
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the gradient vector coefficient at ip.
QuadratureSpaceBase * GetSpace()
Get the associated QuadratureSpaceBase object.
Definition: qfunction.hpp:81
Class for integration point with weight.
Definition: intrules.hpp:25
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:6291
void SetTime(double t)
Set the time for internally stored coefficients.
virtual void SetTime(double t)
Set the time for time dependent coefficients.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the scalar divergence coefficient at ip.
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:348
void ProjectTranspose(MatrixCoefficient &coeff)
Project the transpose of coeff.
Base class for symmetric matrix coefficients that optionally depend on time and space.
DivergenceGridFunctionCoefficient(const GridFunction *gf)
Construct the coefficient with a vector grid function gf. The grid function is not owned by the coeff...
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
void SetTime(double t)
Set the time for internally stored coefficients.
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient.
int dim
Definition: ex24.cpp:53
void SetDirection(const Vector &d_)
QuadratureFunction * qf
Internal QuadratureFunction (owned, may be NULL).
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
Extension of GetGradient(...) for a collection of IntegrationPoints.
Definition: gridfunc.cpp:1715
void SetGridFunction(const GridFunction *gf)
Set the grid function for this coefficient. Also sets the Vector dimension to match that of the gf...
IsoparametricTransformation Transf
Definition: eltrans.hpp:464
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:46
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:804
void SetTime(double t)
Set the time for internally stored coefficients.
MatrixSumCoefficient(MatrixCoefficient &A, MatrixCoefficient &B, double alpha_=1.0, double beta_=1.0)
Construct with the two coefficients. Result is alpha_ * A + beta_ * B.
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:717
virtual void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
Definition: coefficient.cpp:51
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
VectorArrayCoefficient(int dim)
Construct vector of dim coefficients. The actual coefficients still need to be added with Set()...
virtual void Eval(DenseSymmetricMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
DenseSymmetricMatrix mat
Internal matrix used when evaluating this coefficient as a DenseMatrix.
const double alpha
Definition: ex15.cpp:369
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:120
Vector data type.
Definition: vector.hpp:60
virtual void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points.
virtual void SetTime(double t)
Set the time for time dependent coefficients.
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:576
virtual void Project(QuadratureFunction &qf, bool transpose=false)
Fill the QuadratureFunction qf by evaluating the coefficient at the quadrature points. The matrix will be transposed or not according to the boolean argument transpose.
void Project(QuadratureFunction &qf)
Fill the QuadratureFunction qf with the constant value.
Definition: coefficient.cpp:71
CoefficientVector(QuadratureSpaceBase &qs_, CoefficientStorage storage_=CoefficientStorage::FULL)
Create an empty CoefficientVector.
void SetTime(double t)
Set the time for internally stored coefficients.
VectorRotProductCoefficient(VectorCoefficient &A, VectorCoefficient &B)
Constructor with two vector coefficients. Result is .
Defines the coarse-fine transformations of all fine elements.
Definition: ncmesh.hpp:70
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
QuadratureSpaceBase & qs
Associated QuadratureSpaceBase.
void UseExternalData(double *d, int h, int w)
Change the data array and the size of the DenseMatrix.
Definition: densemat.hpp:80
int GetSize() const
Return the total number of quadrature points.
Definition: qspace.hpp:54
virtual void EvalSymmetric(Vector &K, ElementTransformation &T, const IntegrationPoint &ip)
(DEPRECATED) Evaluate the symmetric matrix coefficient at ip.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
void SetTime(double t)
Set the time for internally stored coefficients.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient.
void SetTime(double t)
Set the time for internally stored coefficients.
Class for parallel meshes.
Definition: pmesh.hpp:32
InverseMatrixCoefficient(MatrixCoefficient &A)
Construct with the matrix coefficient. Result is .
void SetTime(double t)
Set the time for internally stored coefficients.
Represents values or vectors of values at quadrature points on a mesh.
Definition: qfunction.hpp:23
void SetTime(double t)
Set the time for internally stored coefficients.
void SetTime(double t)
Set the time for internally stored coefficients.
void SetTime(double t)
Set the time for internally stored coefficients.
ScalarVectorProductCoefficient(double A, VectorCoefficient &B)
Constructor with constant and vector coefficient. Result is A * B.
int GetOrder() const
Return the order of the quadrature rule(s) used by all elements.
Definition: qspace.hpp:57
VectorGridFunctionCoefficient()
Construct an empty coefficient. Calling Eval() before the grid function is set will cause a segfault...
void SetTime(double t)
Set the time for internally stored coefficients.