MFEM  v3.1 Finite element discretization library
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
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  const IntegrationRule &ir)
71 {
72  Vector Mi;
73  M.SetSize(vdim, ir.GetNPoints());
74  for (int i = 0; i < ir.GetNPoints(); i++)
75  {
76  M.GetColumnReference(i, Mi);
77  const IntegrationPoint &ip = ir.IntPoint(i);
78  T.SetIntPoint(&ip);
79  Eval(Mi, T, ip);
80  }
81 }
82
84  const IntegrationPoint &ip)
85 {
86  double x[3];
87  Vector transip(x, 3);
88
89  T.Transform(ip, transip);
90
91  V.SetSize(vdim);
92  if (Function)
93  {
94  (*Function)(transip, V);
95  }
96  else
97  {
98  (*TDFunction)(transip, GetTime(), V);
99  }
100  if (Q)
101  {
102  V *= Q->Eval(T, ip, GetTime());
103  }
104 }
105
107  : VectorCoefficient(dim), Coeff(dim)
108 {
109  for (int i = 0; i < dim; i++)
110  {
111  Coeff[i] = NULL;
112  }
113 }
114
116 {
117  for (int i = 0; i < vdim; i++)
118  {
119  delete Coeff[i];
120  }
121 }
122
124  const IntegrationPoint &ip)
125 {
126  V.SetSize(vdim);
127  for (int i = 0; i < vdim; i++)
128  {
129  V(i) = Coeff[i]->Eval(T, ip, GetTime());
130  }
131 }
132
134  GridFunction *gf) : VectorCoefficient (gf -> VectorDim())
135 {
136  GridFunc = gf;
137 }
138
140  const IntegrationPoint &ip)
141 {
142  GridFunc->GetVectorValue(T.ElementNo, ip, V);
143 }
144
147 {
148  GridFunc->GetVectorValues(T, ir, M);
149 }
150
152  const IntegrationPoint &ip)
153 {
154  V.SetSize(vdim);
155  if (active_attr[T.Attribute-1])
156  {
157  c->SetTime(GetTime());
158  c->Eval(V, T, ip);
159  }
160  else
161  {
162  V = 0.0;
163  }
164 }
165
168 {
169  if (active_attr[T.Attribute-1])
170  {
171  c->SetTime(GetTime());
172  c->Eval(M, T, ir);
173  }
174  else
175  {
176  M.SetSize(vdim);
177  M = 0.0;
178  }
179 }
180
182  const IntegrationPoint &ip)
183 {
184  double x[3];
185  Vector transip(x, 3);
186
187  T.Transform(ip, transip);
188
189  K.SetSize(vdim);
190
191  if (Function)
192  {
193  (*Function)(transip, K);
194  }
195  else
196  {
197  (*TDFunction)(transip, GetTime(), K);
198  }
199 }
200
202  : MatrixCoefficient (dim)
203 {
204  Coeff.SetSize (vdim*vdim);
205 }
206
208 {
209  for (int i=0; i< vdim*vdim; i++)
210  {
211  delete Coeff[i];
212  }
213 }
214
216  const IntegrationPoint &ip)
217 {
218  int i, j;
219
220  for (i = 0; i < vdim; i++)
221  for (j = 0; j < vdim; j++)
222  {
223  K(i,j) = Coeff[i*vdim+j] -> Eval(T, ip, GetTime());
224  }
225 }
226
227 double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh,
228  const IntegrationRule *irs[])
229 {
230  double norm = 0.0;
232
233  for (int i = 0; i < mesh.GetNE(); i++)
234  {
235  tr = mesh.GetElementTransformation(i);
236  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
237  for (int j = 0; j < ir.GetNPoints(); j++)
238  {
239  const IntegrationPoint &ip = ir.IntPoint(j);
240  tr->SetIntPoint(&ip);
241  double val = fabs(coeff.Eval(*tr, ip));
242  if (p < numeric_limits<double>::infinity())
243  {
244  norm += ip.weight * tr->Weight() * pow(val, p);
245  }
246  else
247  {
248  if (norm < val)
249  {
250  norm = val;
251  }
252  }
253  }
254  }
255  return norm;
256 }
257
258 double LpNormLoop(double p, VectorCoefficient &coeff, Mesh &mesh,
259  const IntegrationRule *irs[])
260 {
261  double norm = 0.0;
263  int vdim = coeff.GetVDim();
264  Vector vval(vdim);
265  double val;
266
267  for (int i = 0; i < mesh.GetNE(); i++)
268  {
269  tr = mesh.GetElementTransformation(i);
270  const IntegrationRule &ir = *irs[mesh.GetElementType(i)];
271  for (int j = 0; j < ir.GetNPoints(); j++)
272  {
273  const IntegrationPoint &ip = ir.IntPoint(j);
274  tr->SetIntPoint(&ip);
275  coeff.Eval(vval, *tr, ip);
276  if (p < numeric_limits<double>::infinity())
277  {
278  for (int idim(0); idim < vdim; ++idim)
279  {
280  norm += ip.weight * tr->Weight() * pow(fabs( vval(idim) ), p);
281  }
282  }
283  else
284  {
285  for (int idim(0); idim < vdim; ++idim)
286  {
287  val = fabs(vval(idim));
288  if (norm < val)
289  {
290  norm = val;
291  }
292  }
293  }
294  }
295  }
296
297  return norm;
298 }
299
300 double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh,
301  const IntegrationRule *irs[])
302 {
303  double norm = LpNormLoop(p, coeff, mesh, irs);
304
305  if (p < numeric_limits<double>::infinity())
306  {
307  // negative quadrature weights may cause norm to be negative
308  if (norm < 0.0)
309  {
310  norm = -pow(-norm, 1.0/p);
311  }
312  else
313  {
314  norm = pow(norm, 1.0/p);
315  }
316  }
317
318  return norm;
319 }
320
321 double ComputeLpNorm(double p, VectorCoefficient &coeff, Mesh &mesh,
322  const IntegrationRule *irs[])
323 {
324  double norm = LpNormLoop(p, coeff, mesh, irs);
325
326  if (p < numeric_limits<double>::infinity())
327  {
328  // negative quadrature weights may cause norm to be negative
329  if (norm < 0.0)
330  {
331  norm = -pow(-norm, 1.0/p);
332  }
333  else
334  {
335  norm = pow(norm, 1.0/p);
336  }
337  }
338
339  return norm;
340 }
341
342 #ifdef MFEM_USE_MPI
343 double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh,
344  const IntegrationRule *irs[])
345 {
346  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
347  double glob_norm = 0;
348
349  MPI_Comm comm = pmesh.GetComm();
350
351  if (p < numeric_limits<double>::infinity())
352  {
353  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
354
355  // negative quadrature weights may cause norm to be negative
356  if (glob_norm < 0.0)
357  {
358  glob_norm = -pow(-glob_norm, 1.0/p);
359  }
360  else
361  {
362  glob_norm = pow(glob_norm, 1.0/p);
363  }
364  }
365  else
366  {
367  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
368  }
369
370  return glob_norm;
371 }
372
373 double ComputeGlobalLpNorm(double p, VectorCoefficient &coeff, ParMesh &pmesh,
374  const IntegrationRule *irs[])
375 {
376  double loc_norm = LpNormLoop(p, coeff, pmesh, irs);
377  double glob_norm = 0;
378
379  MPI_Comm comm = pmesh.GetComm();
380
381  if (p < numeric_limits<double>::infinity())
382  {
383  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
384
385  // negative quadrature weights may cause norm to be negative
386  if (glob_norm < 0.0)
387  {
388  glob_norm = -pow(-glob_norm, 1.0/p);
389  }
390  else
391  {
392  glob_norm = pow(glob_norm, 1.0/p);
393  }
394  }
395  else
396  {
397  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
398  }
399
400  return glob_norm;
401 }
402 #endif
403
404 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:228
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Class for integration rule.
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:292
void SetSize(int s)
Resizes the vector if the new size is different.
Definition: vector.hpp:259
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:33
Data type dense matrix using column-major storage.
Definition: densemat.hpp:22
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Definition: coefficient.cpp:83
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:454
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Definition: coefficient.cpp:55
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate coefficient.
Definition: coefficient.cpp:31
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:231
int dim
Definition: ex3.cpp:48
double ComputeLpNorm(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh, const IntegrationRule *irs[])
int GetElementType(int i) const
Returns the type of element i.
Definition: mesh.cpp:4382
double Eval(int i, int j, ElementTransformation &T, IntegrationPoint &ip)
virtual ~VectorArrayCoefficient()
Destroys vector coefficient.
virtual double Weight()=0
int GetVDim()
Returns dimension of the vector.
double Eval(int i, ElementTransformation &T, IntegrationPoint &ip)
Evaluates i&#39;th component of the vector.
MPI_Comm GetComm() const
Definition: pmesh.hpp:86
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:200
void SetTime(double t)
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
double ComputeGlobalLpNorm(double p, Coefficient &coeff, ParMesh &pmesh, const IntegrationRule *irs[])
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:173
void GetVectorValues(ElementTransformation &T, const IntegrationRule &ir, DenseMatrix &vals) const
Definition: gridfunc.cpp:405
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
VectorArrayCoefficient(int dim)
Construct vector of dim coefficients.
Vector data type.
Definition: vector.hpp:33
virtual void Transform(const IntegrationPoint &, Vector &)=0
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:71
Class for parallel meshes.
Definition: pmesh.hpp:28