MFEM  v3.4
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, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 // Implementation of 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.ElementNo, 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 
90 
92  const IntegrationRule &ir)
93 {
94  Vector Mi;
95  M.SetSize(vdim, ir.GetNPoints());
96  for (int i = 0; i < ir.GetNPoints(); i++)
97  {
98  M.GetColumnReference(i, Mi);
99  const IntegrationPoint &ip = ir.IntPoint(i);
100  T.SetIntPoint(&ip);
101  Eval(Mi, T, ip);
102  }
103 }
104 
106  const IntegrationPoint &ip)
107 {
108  double x[3];
109  Vector transip(x, 3);
110 
111  T.Transform(ip, transip);
112 
113  V.SetSize(vdim);
114  if (Function)
115  {
116  (*Function)(transip, V);
117  }
118  else
119  {
120  (*TDFunction)(transip, GetTime(), V);
121  }
122  if (Q)
123  {
124  V *= Q->Eval(T, ip, GetTime());
125  }
126 }
127 
129  : VectorCoefficient(dim), Coeff(dim)
130 {
131  for (int i = 0; i < dim; i++)
132  {
133  Coeff[i] = NULL;
134  }
135 }
136 
138 {
139  for (int i = 0; i < vdim; i++)
140  {
141  delete Coeff[i];
142  }
143 }
144 
146  const IntegrationPoint &ip)
147 {
148  V.SetSize(vdim);
149  for (int i = 0; i < vdim; i++)
150  {
151  V(i) = this->Eval(i, T, ip);
152  }
153 }
154 
156  GridFunction *gf) : VectorCoefficient (gf -> VectorDim())
157 {
158  GridFunc = gf;
159 }
160 
162  const IntegrationPoint &ip)
163 {
164  GridFunc->GetVectorValue(T.ElementNo, ip, V);
165 }
166 
169 {
170  GridFunc->GetVectorValues(T, ir, M);
171 }
172 
174 {
175  dir = _d;
176  (*this).vdim = dir.Size();
177 }
178 
181 {
182  V = dir;
183  d.SetTime(GetTime());
184  V *= d.EvalDelta(T, ip);
185 }
186 
188  const IntegrationPoint &ip)
189 {
190  V.SetSize(vdim);
191  if (active_attr[T.Attribute-1])
192  {
193  c->SetTime(GetTime());
194  c->Eval(V, T, ip);
195  }
196  else
197  {
198  V = 0.0;
199  }
200 }
201 
204 {
205  if (active_attr[T.Attribute-1])
206  {
207  c->SetTime(GetTime());
208  c->Eval(M, T, ir);
209  }
210  else
211  {
212  M.SetSize(vdim);
213  M = 0.0;
214  }
215 }
216 
218  const IntegrationPoint &ip)
219 {
220  double x[3];
221  Vector transip(x, 3);
222 
223  T.Transform(ip, transip);
224 
225  K.SetSize(height, width);
226 
227  if (Function)
228  {
229  (*Function)(transip, K);
230  }
231  else if (TDFunction)
232  {
233  (*TDFunction)(transip, GetTime(), K);
234  }
235  else
236  {
237  K = mat;
238  }
239  if (Q)
240  {
241  K *= Q->Eval(T, ip, GetTime());
242  }
243 }
244 
246  : MatrixCoefficient (dim)
247 {
248  Coeff.SetSize(height*width);
249  for (int i = 0; i < (height*width); i++)
250  {
251  Coeff[i] = NULL;
252  }
253 }
254 
256 {
257  for (int i=0; i < height*width; i++)
258  {
259  delete Coeff[i];
260  }
261 }
262 
264  const IntegrationPoint &ip)
265 {
266  for (int i = 0; i < height; i++)
267  {
268  for (int j = 0; j < width; j++)
269  {
270  K(i,j) = this->Eval(i, j, T, ip);
271  }
272  }
273 }
274 
276  const IntegrationPoint &ip)
277 {
278  if (active_attr[T.Attribute-1])
279  {
280  c->SetTime(GetTime());
281  c->Eval(K, T, ip);
282  }
283  else
284  {
285  K.SetSize(height, width);
286  K = 0.0;
287  }
288 }
289 
290 double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh,
291  const IntegrationRule *irs[])
292 {
293  double norm = 0.0;
295 
296  for (int i = 0; i < mesh.GetNE(); i++)
297  {
298  tr = mesh.GetElementTransformation(i);
299  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
300  for (int j = 0; j < ir.GetNPoints(); j++)
301  {
302  const IntegrationPoint &ip = ir.IntPoint(j);
303  tr->SetIntPoint(&ip);
304  double val = fabs(coeff.Eval(*tr, ip));
305  if (p < infinity())
306  {
307  norm += ip.weight * tr->Weight() * pow(val, p);
308  }
309  else
310  {
311  if (norm < val)
312  {
313  norm = val;
314  }
315  }
316  }
317  }
318  return norm;
319 }
320 
321 double LpNormLoop(double p, VectorCoefficient &coeff, Mesh &mesh,
322  const IntegrationRule *irs[])
323 {
324  double norm = 0.0;
326  int vdim = coeff.GetVDim();
327  Vector vval(vdim);
328  double val;
329 
330  for (int i = 0; i < mesh.GetNE(); i++)
331  {
332  tr = mesh.GetElementTransformation(i);
333  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
334  for (int j = 0; j < ir.GetNPoints(); j++)
335  {
336  const IntegrationPoint &ip = ir.IntPoint(j);
337  tr->SetIntPoint(&ip);
338  coeff.Eval(vval, *tr, ip);
339  if (p < infinity())
340  {
341  for (int idim(0); idim < vdim; ++idim)
342  {
343  norm += ip.weight * tr->Weight() * pow(fabs( vval(idim) ), p);
344  }
345  }
346  else
347  {
348  for (int idim(0); idim < vdim; ++idim)
349  {
350  val = fabs(vval(idim));
351  if (norm < val)
352  {
353  norm = val;
354  }
355  }
356  }
357  }
358  }
359 
360  return norm;
361 }
362 
363 double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh,
364  const IntegrationRule *irs[])
365 {
366  double norm = LpNormLoop(p, coeff, mesh, irs);
367 
368  if (p < infinity())
369  {
370  // negative quadrature weights may cause norm to be negative
371  if (norm < 0.0)
372  {
373  norm = -pow(-norm, 1.0/p);
374  }
375  else
376  {
377  norm = pow(norm, 1.0/p);
378  }
379  }
380 
381  return norm;
382 }
383 
384 double ComputeLpNorm(double p, VectorCoefficient &coeff, Mesh &mesh,
385  const IntegrationRule *irs[])
386 {
387  double norm = LpNormLoop(p, coeff, mesh, irs);
388 
389  if (p < infinity())
390  {
391  // negative quadrature weights may cause norm to be negative
392  if (norm < 0.0)
393  {
394  norm = -pow(-norm, 1.0/p);
395  }
396  else
397  {
398  norm = pow(norm, 1.0/p);
399  }
400  }
401 
402  return norm;
403 }
404 
405 #ifdef MFEM_USE_MPI
406 double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh,
407  const IntegrationRule *irs[])
408 {
409  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
410  double glob_norm = 0;
411 
412  MPI_Comm comm = pmesh.GetComm();
413 
414  if (p < infinity())
415  {
416  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
417 
418  // negative quadrature weights may cause norm to be negative
419  if (glob_norm < 0.0)
420  {
421  glob_norm = -pow(-glob_norm, 1.0/p);
422  }
423  else
424  {
425  glob_norm = pow(glob_norm, 1.0/p);
426  }
427  }
428  else
429  {
430  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
431  }
432 
433  return glob_norm;
434 }
435 
436 double ComputeGlobalLpNorm(double p, VectorCoefficient &coeff, ParMesh &pmesh,
437  const IntegrationRule *irs[])
438 {
439  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
440  double glob_norm = 0;
441 
442  MPI_Comm comm = pmesh.GetComm();
443 
444  if (p < infinity())
445  {
446  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
447 
448  // negative quadrature weights may cause norm to be negative
449  if (glob_norm < 0.0)
450  {
451  glob_norm = -pow(-glob_norm, 1.0/p);
452  }
453  else
454  {
455  glob_norm = pow(glob_norm, 1.0/p);
456  }
457  }
458  else
459  {
460  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
461  }
462 
463  return glob_norm;
464 }
465 #endif
466 
467 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:232
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
VectorGridFunctionCoefficient(GridFunction *gf)
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:400
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:328
double Eval(int i, ElementTransformation &T, const IntegrationPoint &ip)
Evaluates i&#39;th component of the vector.
void SetDeltaCenter(const Vector &center)
Definition: coefficient.cpp:69
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:120
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:618
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Definition: coefficient.cpp:55
void SetTime(double t)
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate coefficient.
Definition: coefficient.cpp:31
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
double Eval(int i, int j, ElementTransformation &T, const IntegrationPoint &ip)
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:235
int dim
Definition: ex3.cpp:47
double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
void SetDirection(const Vector &_d)
int GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:4029
virtual ~VectorArrayCoefficient()
Destroys vector coefficient.
int GetVDim()
Returns dimension of the vector.
void SetTime(double t)
Definition: coefficient.hpp:39
void GetDeltaCenter(Vector &center)
Definition: coefficient.cpp:77
MPI_Comm GetComm() const
Definition: pmesh.hpp:123
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient function.
Definition: coefficient.cpp:24
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:243
void SetTime(double t)
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
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[])
double EvalDelta(ElementTransformation &T, const IntegrationPoint &ip)
Return the Scale() multiplied by the weight Coefficient, if any.
Definition: coefficient.cpp:83
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Definition: coefficient.cpp:49
Class for integration point with weight.
Definition: intrules.hpp:25
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:279
void GetVectorValues(ElementTransformation &T, const IntegrationRule &ir, DenseMatrix &vals) const
Definition: gridfunc.cpp:517
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:42
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
VectorArrayCoefficient(int dim)
Construct vector of dim coefficients.
Vector data type.
Definition: vector.hpp:48
virtual void Transform(const IntegrationPoint &, Vector &)=0
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:86
Class for parallel meshes.
Definition: pmesh.hpp:32