MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
coefficient.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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 
25  const IntegrationPoint & ip)
26 {
27  int att = T.Attribute;
28  return (constants(att-1));
29 }
30 
32  const IntegrationPoint & ip)
33 {
34  double x[3];
35  Vector transip(x, 3);
36 
37  T.Transform(ip, transip);
38 
39  if (Function)
40  {
41  return Function(transip);
42  }
43  else
44  {
45  return TDFunction(transip, GetTime());
46  }
47 }
48 
50  const IntegrationPoint &ip)
51 {
52  return GridF -> GetValue (T, ip, Component);
53 }
54 
56  const IntegrationPoint &ip)
57 {
58  if (Q2)
59  {
60  return (*Transform2)(Q1->Eval(T, ip, GetTime()),
61  Q2->Eval(T, ip, GetTime()));
62  }
63  else
64  {
65  return (*Transform1)(Q1->Eval(T, ip, GetTime()));
66  }
67 }
68 
70 {
71  MFEM_VERIFY(vcenter.Size() <= 3,
72  "SetDeltaCenter::Maximum number of dim supported is 3")
73  for (int i = 0; i < vcenter.Size(); i++) { center[i] = vcenter[i]; }
74  sdim = vcenter.Size();
75 }
76 
78 {
79  vcenter.SetSize(sdim);
80  vcenter = center;
81 }
82 
84  const IntegrationPoint &ip)
85 {
86  double w = Scale();
87  return weight ? weight->Eval(T, ip, GetTime())*w : w;
88 }
89 
91  const IntegrationRule &ir)
92 {
93  Vector Mi;
94  M.SetSize(vdim, ir.GetNPoints());
95  for (int i = 0; i < ir.GetNPoints(); i++)
96  {
97  M.GetColumnReference(i, Mi);
98  const IntegrationPoint &ip = ir.IntPoint(i);
99  T.SetIntPoint(&ip);
100  Eval(Mi, T, ip);
101  }
102 }
103 
105  const IntegrationPoint &ip)
106 {
107  double x[3];
108  Vector transip(x, 3);
109 
110  T.Transform(ip, transip);
111 
112  V.SetSize(vdim);
113  if (Function)
114  {
115  Function(transip, V);
116  }
117  else
118  {
119  TDFunction(transip, GetTime(), V);
120  }
121  if (Q)
122  {
123  V *= Q->Eval(T, ip, GetTime());
124  }
125 }
126 
128  : VectorCoefficient(dim), Coeff(dim), ownCoeff(dim)
129 {
130  for (int i = 0; i < dim; i++)
131  {
132  Coeff[i] = NULL;
133  ownCoeff[i] = true;
134  }
135 }
136 
137 void VectorArrayCoefficient::Set(int i, Coefficient *c, bool own)
138 {
139  if (ownCoeff[i]) { delete Coeff[i]; }
140  Coeff[i] = c;
141  ownCoeff[i] = own;
142 }
143 
145 {
146  for (int i = 0; i < vdim; i++)
147  {
148  if (ownCoeff[i]) { delete Coeff[i]; }
149  }
150 }
151 
153  const IntegrationPoint &ip)
154 {
155  V.SetSize(vdim);
156  for (int i = 0; i < vdim; i++)
157  {
158  V(i) = this->Eval(i, T, ip);
159  }
160 }
161 
163  const GridFunction *gf)
164  : VectorCoefficient ((gf) ? gf -> VectorDim() : 0)
165 {
166  GridFunc = gf;
167 }
168 
170 {
171  GridFunc = gf; vdim = (gf) ? gf -> VectorDim() : 0;
172 }
173 
175  const IntegrationPoint &ip)
176 {
177  GridFunc->GetVectorValue(T, ip, V);
178 }
179 
182 {
183  GridFunc->GetVectorValues(T, ir, M);
184 }
185 
187  const GridFunction *gf)
188  : VectorCoefficient((gf) ?
189  gf -> FESpace() -> GetMesh() -> SpaceDimension() : 0)
190 {
191  GridFunc = gf;
192 }
193 
195 {
196  GridFunc = gf; vdim = (gf) ?
197  gf -> FESpace() -> GetMesh() -> SpaceDimension() : 0;
198 }
199 
201  const IntegrationPoint &ip)
202 {
203  GridFunc->GetGradient(T, V);
204 }
205 
208 {
209  GridFunc->GetGradients(T, ir, M);
210 }
211 
213  const GridFunction *gf)
214  : VectorCoefficient(0)
215 {
216  SetGridFunction(gf);
217 }
218 
220 {
221  if (gf)
222  {
223  int sdim = gf -> FESpace() -> GetMesh() -> SpaceDimension();
224  MFEM_VERIFY(sdim == 2 || sdim == 3,
225  "CurlGridFunctionCoefficient "
226  "only defind for spaces of dimension 2 or 3.");
227  }
228  GridFunc = gf;
229  vdim = (gf) ? (2 * gf -> FESpace() -> GetMesh() -> SpaceDimension() - 3) : 0;
230 }
231 
233  const IntegrationPoint &ip)
234 {
235  GridFunc->GetCurl(T, V);
236 }
237 
239  const GridFunction *gf) : Coefficient()
240 {
241  GridFunc = gf;
242 }
243 
245  const IntegrationPoint &ip)
246 {
247  return GridFunc->GetDivergence(T);
248 }
249 
251 {
252  dir = _d;
253  (*this).vdim = dir.Size();
254 }
255 
258 {
259  V = dir;
260  d.SetTime(GetTime());
261  V *= d.EvalDelta(T, ip);
262 }
263 
265  const IntegrationPoint &ip)
266 {
267  V.SetSize(vdim);
268  if (active_attr[T.Attribute-1])
269  {
270  c->SetTime(GetTime());
271  c->Eval(V, T, ip);
272  }
273  else
274  {
275  V = 0.0;
276  }
277 }
278 
281 {
282  if (active_attr[T.Attribute-1])
283  {
284  c->SetTime(GetTime());
285  c->Eval(M, T, ir);
286  }
287  else
288  {
289  M.SetSize(vdim, ir.GetNPoints());
290  M = 0.0;
291  }
292 }
293 
295  const IntegrationPoint &ip)
296 {
297  double x[3];
298  Vector transip(x, 3);
299 
300  T.Transform(ip, transip);
301 
302  K.SetSize(height, width);
303 
304  if (symmetric) // Use SymmFunction
305  {
306  MFEM_VERIFY(height == width && SymmFunction,
307  "MatrixFunctionCoefficient is not symmetric");
308 
309  Vector Ksym((width * (width + 1)) / 2); // 1x1: 1, 2x2: 3, 3x3: 6
310 
311  SymmFunction(transip, Ksym);
312 
313  // Copy upper triangular values from Ksym to the full matrix K
314  int os = 0;
315  for (int i=0; i<height; ++i)
316  {
317  for (int j=i; j<width; ++j)
318  {
319  const double Kij = Ksym[j - i + os];
320  K(i,j) = Kij;
321  if (j != i) { K(j,i) = Kij; }
322  }
323 
324  os += width - i;
325  }
326  }
327  else
328  {
329  if (Function)
330  {
331  Function(transip, K);
332  }
333  else if (TDFunction)
334  {
335  TDFunction(transip, GetTime(), K);
336  }
337  else
338  {
339  K = mat;
340  }
341  }
342 
343  if (Q)
344  {
345  K *= Q->Eval(T, ip, GetTime());
346  }
347 }
348 
351  const IntegrationPoint &ip)
352 {
353  MFEM_VERIFY(symmetric && height == width && SymmFunction,
354  "MatrixFunctionCoefficient is not symmetric");
355 
356  double x[3];
357  Vector transip(x, 3);
358 
359  T.Transform(ip, transip);
360 
361  K.SetSize((width * (width + 1)) / 2); // 1x1: 1, 2x2: 3, 3x3: 6
362 
363  if (SymmFunction)
364  {
365  SymmFunction(transip, K);
366  }
367 
368  if (Q)
369  {
370  K *= Q->Eval(T, ip, GetTime());
371  }
372 }
373 
375  : MatrixCoefficient (dim)
376 {
377  Coeff.SetSize(height*width);
378  ownCoeff.SetSize(height*width);
379  for (int i = 0; i < (height*width); i++)
380  {
381  Coeff[i] = NULL;
382  ownCoeff[i] = true;
383  }
384 }
385 
386 void MatrixArrayCoefficient::Set(int i, int j, Coefficient * c, bool own)
387 {
388  if (ownCoeff[i*width+j]) { delete Coeff[i*width+j]; }
389  Coeff[i*width+j] = c;
390  ownCoeff[i*width+j] = own;
391 }
392 
394 {
395  for (int i=0; i < height*width; i++)
396  {
397  if (ownCoeff[i]) { delete Coeff[i]; }
398  }
399 }
400 
402  const IntegrationPoint &ip)
403 {
404  for (int i = 0; i < height; i++)
405  {
406  for (int j = 0; j < width; j++)
407  {
408  K(i,j) = this->Eval(i, j, T, ip);
409  }
410  }
411 }
412 
414  const IntegrationPoint &ip)
415 {
416  if (active_attr[T.Attribute-1])
417  {
418  c->SetTime(GetTime());
419  c->Eval(K, T, ip);
420  }
421  else
422  {
423  K.SetSize(height, width);
424  K = 0.0;
425  }
426 }
427 
430  : a(&A), b(&B)
431 {
432  MFEM_ASSERT(A.GetVDim() == B.GetVDim(),
433  "InnerProductCoefficient: "
434  "Arguments have incompatible dimensions.");
435 }
436 
438  const IntegrationPoint &ip)
439 {
440  a->Eval(va, T, ip);
441  b->Eval(vb, T, ip);
442  return va * vb;
443 }
444 
447  : a(&A), b(&B), va(A.GetVDim()), vb(B.GetVDim())
448 {
449  MFEM_ASSERT(A.GetVDim() == 2 && B.GetVDim() == 2,
450  "VectorRotProductCoefficient: "
451  "Arguments must have dimension equal to two.");
452 }
453 
455  const IntegrationPoint &ip)
456 {
457  a->Eval(va, T, ip);
458  b->Eval(vb, T, ip);
459  return va[0] * vb[1] - va[1] * vb[0];
460 }
461 
463  : a(&A), ma(A.GetHeight(), A.GetWidth())
464 {
465  MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
466  "DeterminantCoefficient: "
467  "Argument must be a square matrix.");
468 }
469 
471  const IntegrationPoint &ip)
472 {
473  a->Eval(ma, T, ip);
474  return ma.Det();
475 }
476 
478  : VectorCoefficient(dim),
479  ACoef(NULL), BCoef(NULL),
480  A(dim), B(dim),
481  alphaCoef(NULL), betaCoef(NULL),
482  alpha(1.0), beta(1.0)
483 {
484  A = 0.0; B = 0.0;
485 }
486 
488  VectorCoefficient &_B,
489  double _alpha, double _beta)
490  : VectorCoefficient(_A.GetVDim()),
491  ACoef(&_A), BCoef(&_B),
492  A(_A.GetVDim()), B(_A.GetVDim()),
493  alphaCoef(NULL), betaCoef(NULL),
494  alpha(_alpha), beta(_beta)
495 {
496  MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
497  "VectorSumCoefficient: "
498  "Arguments must have the same dimension.");
499 }
500 
502  VectorCoefficient &_B,
503  Coefficient &_alpha,
504  Coefficient &_beta)
505  : VectorCoefficient(_A.GetVDim()),
506  ACoef(&_A), BCoef(&_B),
507  A(_A.GetVDim()),
508  B(_A.GetVDim()),
509  alphaCoef(&_alpha),
510  betaCoef(&_beta),
511  alpha(0.0), beta(0.0)
512 {
513  MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
514  "VectorSumCoefficient: "
515  "Arguments must have the same dimension.");
516 }
517 
519  const IntegrationPoint &ip)
520 {
521  V.SetSize(A.Size());
522  if ( ACoef) { ACoef->Eval(A, T, ip); }
523  if ( BCoef) { BCoef->Eval(B, T, ip); }
524  if (alphaCoef) { alpha = alphaCoef->Eval(T, ip); }
525  if ( betaCoef) { beta = betaCoef->Eval(T, ip); }
526  add(alpha, A, beta, B, V);
527 }
528 
530  double A,
532  : VectorCoefficient(B.GetVDim()), aConst(A), a(NULL), b(&B)
533 {}
534 
536  Coefficient &A,
538  : VectorCoefficient(B.GetVDim()), aConst(0.0), a(&A), b(&B)
539 {}
540 
542  const IntegrationPoint &ip)
543 {
544  double sa = (a == NULL) ? aConst : a->Eval(T, ip);
545  b->Eval(V, T, ip);
546  V *= sa;
547 }
548 
550  double _tol)
551  : VectorCoefficient(A.GetVDim()), a(&A), tol(_tol)
552 {}
553 
555  const IntegrationPoint &ip)
556 {
557  a->Eval(V, T, ip);
558  double nv = V.Norml2();
559  V *= (nv > tol) ? (1.0/nv) : 0.0;
560 }
561 
565  : VectorCoefficient(3), a(&A), b(&B), va(A.GetVDim()), vb(B.GetVDim())
566 {
567  MFEM_ASSERT(A.GetVDim() == 3 && B.GetVDim() == 3,
568  "VectorCrossProductCoefficient: "
569  "Arguments must have dimension equal to three.");
570 }
571 
573  const IntegrationPoint &ip)
574 {
575  a->Eval(va, T, ip);
576  b->Eval(vb, T, ip);
577  V.SetSize(3);
578  V[0] = va[1] * vb[2] - va[2] * vb[1];
579  V[1] = va[2] * vb[0] - va[0] * vb[2];
580  V[2] = va[0] * vb[1] - va[1] * vb[0];
581 }
582 
585  : VectorCoefficient(A.GetHeight()), a(&A), b(&B),
586  ma(A.GetHeight(), A.GetWidth()), vb(B.GetVDim())
587 {
588  MFEM_ASSERT(A.GetWidth() == B.GetVDim(),
589  "MatrixVectorProductCoefficient: "
590  "Arguments have incompatible dimensions.");
591 }
592 
594  const IntegrationPoint &ip)
595 {
596  a->Eval(ma, T, ip);
597  b->Eval(vb, T, ip);
598  ma.Mult(vb, V);
599 }
600 
602  const IntegrationPoint &ip)
603 {
604  M.SetSize(dim);
605  M = 0.0;
606  for (int d=0; d<dim; d++) { M(d,d) = 1.0; }
607 }
608 
611  double _alpha, double _beta)
612  : MatrixCoefficient(A.GetHeight(), A.GetWidth()),
613  a(&A), b(&B), alpha(_alpha), beta(_beta),
614  ma(A.GetHeight(), A.GetWidth())
615 {
616  MFEM_ASSERT(A.GetHeight() == B.GetHeight() && A.GetWidth() == B.GetWidth(),
617  "MatrixSumCoefficient: "
618  "Arguments must have the same dimensions.");
619 }
620 
622  const IntegrationPoint &ip)
623 {
624  b->Eval(M, T, ip);
625  if ( beta != 1.0 ) { M *= beta; }
626  a->Eval(ma, T, ip);
627  M.Add(alpha, ma);
628 }
629 
631  double A,
633  : MatrixCoefficient(B.GetHeight(), B.GetWidth()), aConst(A), a(NULL), b(&B)
634 {}
635 
637  Coefficient &A,
639  : MatrixCoefficient(B.GetHeight(), B.GetWidth()), aConst(0.0), a(&A), b(&B)
640 {}
641 
644  const IntegrationPoint &ip)
645 {
646  double sa = (a == NULL) ? aConst : a->Eval(T, ip);
647  b->Eval(M, T, ip);
648  M *= sa;
649 }
650 
652  : MatrixCoefficient(A.GetWidth(), A.GetHeight()), a(&A)
653 {}
654 
657  const IntegrationPoint &ip)
658 {
659  a->Eval(M, T, ip);
660  M.Transpose();
661 }
662 
664  : MatrixCoefficient(A.GetHeight(), A.GetWidth()), a(&A)
665 {
666  MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
667  "InverseMatrixCoefficient: "
668  "Argument must be a square matrix.");
669 }
670 
673  const IntegrationPoint &ip)
674 {
675  a->Eval(M, T, ip);
676  M.Invert();
677 }
678 
681  : MatrixCoefficient(A.GetVDim(), B.GetVDim()), a(&A), b(&B),
682  va(A.GetVDim()), vb(B.GetVDim())
683 {}
684 
686  const IntegrationPoint &ip)
687 {
688  a->Eval(va, T, ip);
689  b->Eval(vb, T, ip);
690  M.SetSize(va.Size(), vb.Size());
691  for (int i=0; i<va.Size(); i++)
692  {
693  for (int j=0; j<vb.Size(); j++)
694  {
695  M(i, j) = va[i] * vb[j];
696  }
697  }
698 }
699 
702  : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(0.0), a(&A), k(&K),
703  vk(K.GetVDim())
704 {}
705 
707  const IntegrationPoint &ip)
708 {
709  k->Eval(vk, T, ip);
710  M.SetSize(vk.Size(), vk.Size());
711  M = 0.0;
712  double k2 = vk*vk;
713  for (int i=0; i<vk.Size(); i++)
714  {
715  M(i, i) = k2;
716  for (int j=0; j<vk.Size(); j++)
717  {
718  M(i, j) -= vk[i] * vk[j];
719  }
720  }
721  M *= ((a == NULL ) ? aConst : a->Eval(T, ip) );
722 }
723 
724 double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh,
725  const IntegrationRule *irs[])
726 {
727  double norm = 0.0;
729 
730  for (int i = 0; i < mesh.GetNE(); i++)
731  {
732  tr = mesh.GetElementTransformation(i);
733  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
734  for (int j = 0; j < ir.GetNPoints(); j++)
735  {
736  const IntegrationPoint &ip = ir.IntPoint(j);
737  tr->SetIntPoint(&ip);
738  double val = fabs(coeff.Eval(*tr, ip));
739  if (p < infinity())
740  {
741  norm += ip.weight * tr->Weight() * pow(val, p);
742  }
743  else
744  {
745  if (norm < val)
746  {
747  norm = val;
748  }
749  }
750  }
751  }
752  return norm;
753 }
754 
755 double LpNormLoop(double p, VectorCoefficient &coeff, Mesh &mesh,
756  const IntegrationRule *irs[])
757 {
758  double norm = 0.0;
760  int vdim = coeff.GetVDim();
761  Vector vval(vdim);
762  double val;
763 
764  for (int i = 0; i < mesh.GetNE(); i++)
765  {
766  tr = mesh.GetElementTransformation(i);
767  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
768  for (int j = 0; j < ir.GetNPoints(); j++)
769  {
770  const IntegrationPoint &ip = ir.IntPoint(j);
771  tr->SetIntPoint(&ip);
772  coeff.Eval(vval, *tr, ip);
773  if (p < infinity())
774  {
775  for (int idim(0); idim < vdim; ++idim)
776  {
777  norm += ip.weight * tr->Weight() * pow(fabs( vval(idim) ), p);
778  }
779  }
780  else
781  {
782  for (int idim(0); idim < vdim; ++idim)
783  {
784  val = fabs(vval(idim));
785  if (norm < val)
786  {
787  norm = val;
788  }
789  }
790  }
791  }
792  }
793 
794  return norm;
795 }
796 
797 double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh,
798  const IntegrationRule *irs[])
799 {
800  double norm = LpNormLoop(p, coeff, mesh, irs);
801 
802  if (p < infinity())
803  {
804  // negative quadrature weights may cause norm to be negative
805  if (norm < 0.0)
806  {
807  norm = -pow(-norm, 1.0/p);
808  }
809  else
810  {
811  norm = pow(norm, 1.0/p);
812  }
813  }
814 
815  return norm;
816 }
817 
818 double ComputeLpNorm(double p, VectorCoefficient &coeff, Mesh &mesh,
819  const IntegrationRule *irs[])
820 {
821  double norm = LpNormLoop(p, coeff, mesh, irs);
822 
823  if (p < infinity())
824  {
825  // negative quadrature weights may cause norm to be negative
826  if (norm < 0.0)
827  {
828  norm = -pow(-norm, 1.0/p);
829  }
830  else
831  {
832  norm = pow(norm, 1.0/p);
833  }
834  }
835 
836  return norm;
837 }
838 
839 #ifdef MFEM_USE_MPI
840 double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh,
841  const IntegrationRule *irs[])
842 {
843  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
844  double glob_norm = 0;
845 
846  MPI_Comm comm = pmesh.GetComm();
847 
848  if (p < infinity())
849  {
850  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
851 
852  // negative quadrature weights may cause norm to be negative
853  if (glob_norm < 0.0)
854  {
855  glob_norm = -pow(-glob_norm, 1.0/p);
856  }
857  else
858  {
859  glob_norm = pow(glob_norm, 1.0/p);
860  }
861  }
862  else
863  {
864  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
865  }
866 
867  return glob_norm;
868 }
869 
870 double ComputeGlobalLpNorm(double p, VectorCoefficient &coeff, ParMesh &pmesh,
871  const IntegrationRule *irs[])
872 {
873  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
874  double glob_norm = 0;
875 
876  MPI_Comm comm = pmesh.GetComm();
877 
878  if (p < infinity())
879  {
880  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
881 
882  // negative quadrature weights may cause norm to be negative
883  if (glob_norm < 0.0)
884  {
885  glob_norm = -pow(-glob_norm, 1.0/p);
886  }
887  else
888  {
889  glob_norm = pow(glob_norm, 1.0/p);
890  }
891  }
892  else
893  {
894  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
895  }
896 
897  return glob_norm;
898 }
899 #endif
900 
902  QuadratureFunction &qf)
903  : VectorCoefficient(qf.GetVDim()), QuadF(qf), index(0) { }
904 
906 {
907  MFEM_VERIFY(_index >= 0, "Index must be >= 0");
908  MFEM_VERIFY(_index < QuadF.GetVDim(),
909  "Index must be < QuadratureFunction length");
910  index = _index;
911 
912  MFEM_VERIFY(_length > 0, "Length must be > 0");
913  MFEM_VERIFY(_length <= QuadF.GetVDim() - index,
914  "Length must be <= (QuadratureFunction length - index)");
915 
916  vdim = _length;
917 }
918 
921  const IntegrationPoint &ip)
922 {
923  QuadF.HostRead();
924 
925  if (index == 0 && vdim == QuadF.GetVDim())
926  {
927  QuadF.GetElementValues(T.ElementNo, ip.index, V);
928  }
929  else
930  {
931  Vector temp;
932  QuadF.GetElementValues(T.ElementNo, ip.index, temp);
933  V.SetSize(vdim);
934  for (int i = 0; i < vdim; i++)
935  {
936  V(i) = temp(index + i);
937  }
938  }
939 
940  return;
941 }
942 
944  QuadratureFunction &qf) : QuadF(qf)
945 {
946  MFEM_VERIFY(qf.GetVDim() == 1, "QuadratureFunction's vdim must be 1");
947 }
948 
950  const IntegrationPoint &ip)
951 {
952  QuadF.HostRead();
953  Vector temp(1);
954  QuadF.GetElementValues(T.ElementNo, ip.index, temp);
955  return temp[0];
956 }
957 
958 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
int GetVDim() const
Get the vector dimension.
Definition: gridfunc.hpp:738
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient at ip.
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:633
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 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.
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:414
NormalizedVectorCoefficient(VectorCoefficient &A, double tol=1e-6)
Return a vector normalized to a length of one.
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:459
MatrixVectorProductCoefficient(MatrixCoefficient &A, VectorCoefficient &B)
Constructor with two coefficients. Result is A*B.
double Eval(int i, ElementTransformation &T, const IntegrationPoint &ip)
double Det() const
Definition: densemat.cpp:451
void SetDeltaCenter(const Vector &center)
Set the center location of the delta function.
Definition: coefficient.cpp:69
double GetTime()
Get the time for time dependent coefficients.
Element::Type GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:5030
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:744
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 ...
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:160
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...
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:737
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.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Definition: coefficient.cpp:55
void SetTime(double t)
Set the time for time dependent coefficients.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Definition: coefficient.cpp:31
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:261
void SetComponent(int _index, int _length)
DeterminantCoefficient(MatrixCoefficient &A)
Construct with the matrix.
ScalarMatrixProductCoefficient(double A, MatrixCoefficient &B)
Constructor with one coefficient. Result is A*B.
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.
int GetHeight() const
Get the height of the matrix.
double Eval(int i, int j, ElementTransformation &T, const IntegrationPoint &ip)
const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:376
InnerProductCoefficient(VectorCoefficient &A, VectorCoefficient &B)
Construct with the two vector coefficients. Result is .
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:248
double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
Compute the Lp norm of a function f. .
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.
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:544
double b
Definition: lissajous.cpp:42
void SetDirection(const Vector &_d)
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:637
void Set(int i, Coefficient *c, bool own=true)
Sets coefficient in the vector.
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
virtual ~VectorArrayCoefficient()
Destroys vector coefficient.
int GetVDim()
Returns dimension of the vector.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
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 time dependent coefficients.
Definition: coefficient.hpp:48
void GetDeltaCenter(Vector &center)
Write the center of the delta function into center.
Definition: coefficient.cpp:77
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:1378
MPI_Comm GetComm() const
Definition: pmesh.hpp:236
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
VectorQuadratureFunctionCoefficient(QuadratureFunction &qf)
Constructor with a quadrature function as input.
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:654
Base class for Matrix Coefficients that optionally depend on time and space.
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient.
Definition: coefficient.cpp:24
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:270
void SetTime(double t)
Set the time for time dependent coefficients.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient at ip.
virtual void EvalDelta(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Return the specified direction vector multiplied by the value returned by DeltaCoefficient::EvalDelta...
double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
Compute the global Lp norm of a function f. .
double GetDivergence(ElementTransformation &tr) const
Definition: gridfunc.cpp:1359
virtual double EvalDelta(ElementTransformation &T, const IntegrationPoint &ip)
The value of the function assuming we are evaluating at the delta center.
Definition: coefficient.cpp:83
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Definition: coefficient.cpp:49
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.
Class for integration point with weight.
Definition: intrules.hpp:25
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the scalar divergence coefficient at ip.
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:336
MatrixSumCoefficient(MatrixCoefficient &A, MatrixCoefficient &B, double _alpha=1.0, double _beta=1.0)
Construct with the two coefficients. Result is _alpha * A + _beta * B.
DivergenceGridFunctionCoefficient(const GridFunction *gf)
Construct the coefficient with a vector grid function gf. The grid function is not owned by the coeff...
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the matrix coefficient at ip.
int dim
Definition: ex24.cpp:53
void GetElementValues(int idx, Vector &values)
Return all values associated with mesh element idx in a Vector.
Definition: gridfunc.hpp:883
void SetGridFunction(const GridFunction *gf)
Set the grid function for this coefficient. Also sets the Vector dimension to match that of the gf...
int index(int i, int j, int nx, int ny)
Definition: life.cpp:241
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
Definition: gridfunc.cpp:1617
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&lt;double&gt;::infinity()
Definition: vector.hpp:45
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()...
const double alpha
Definition: ex15.cpp:336
void GetCurl(ElementTransformation &tr, Vector &curl) const
Definition: gridfunc.cpp:1445
Vector data type.
Definition: vector.hpp:51
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:175
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
int GetWidth() const
Get the width of the matrix.
VectorRotProductCoefficient(VectorCoefficient &A, VectorCoefficient &B)
Constructor with two vector coefficients. Result is .
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.
virtual void EvalSymmetric(Vector &K, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the symmetric matrix coefficient at ip.
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:90
Class for parallel meshes.
Definition: pmesh.hpp:32
InverseMatrixCoefficient(MatrixCoefficient &A)
Construct with the matrix coefficient. Result is .
Class representing a function through its values (scalar or vector) at quadrature points...
Definition: gridfunc.hpp:670
void GetGradient(ElementTransformation &tr, Vector &grad) const
Definition: gridfunc.cpp:1547
ScalarVectorProductCoefficient(double A, VectorCoefficient &B)
Constructor with constant and vector coefficient. Result is A * B.
VectorGridFunctionCoefficient()
Construct an empty coefficient. Calling Eval() before the grid function is set will cause a segfault...