MFEM  v4.6.0
Finite element discretization library
tcoefficient.hpp
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 #ifndef MFEM_TEMPLATE_COEFFICIENT
13 #define MFEM_TEMPLATE_COEFFICIENT
14 
15 #include "../config/tconfig.hpp"
16 #include "../linalg/ttensor.hpp"
17 #include "../linalg/tlayout.hpp"
18 #include "../linalg/vector.hpp"
19 #include "gridfunc.hpp"
20 
21 namespace mfem
22 {
23 
24 /// Templated coefficient classes, cf. coefficient.?pp
25 
27 {
28 public:
29  static const int rank = 0; // 0 - scalar, 1 - vector, 2 - matrix
30  static const bool is_const = false;
31  static const bool uses_coordinates = false;
32  static const bool uses_Jacobians = false;
33  static const bool uses_attributes = false;
34  static const bool uses_element_idxs = false;
35 };
36 
37 template <typename complex_t = double>
39 {
40 public:
41  static const bool is_const = true;
42  typedef complex_t complex_type;
43 
44  complex_t value;
45 
46  TConstantCoefficient(complex_t val) : value(val) { }
47  // default copy constructor
48 
49  // T_result_t is the transformation result type (not used here).
50  template <typename T_result_t, typename c_layout_t, typename c_data_t>
51  inline MFEM_ALWAYS_INLINE
52  void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c) const
53  {
54  TAssign<AssignOp::Set>(l, c, value);
55  }
56 };
57 
58 
59 /** @brief Function coefficient.
60  @tparam Func has to implement at least one of the following methods,
61  depending on the dimension that will be used:
62  complex_t Eval1D(real_t);
63  complex_t Eval2D(real_t,real_t);
64  complex_t Eval3D(real_t,real_t,real_t);
65  Use MFEM_FLOPS_ADD() to count flops inside Eval*D. */
66 template <typename Func, typename complex_t = double>
68 {
69 public:
70  static const bool uses_coordinates = true;
71  typedef complex_t complex_type;
72 
73 protected:
74  Func F;
75 
76  template <int dim, bool dummy> struct Dim;
77  template <bool dummy> struct Dim<1,dummy>
78  {
79  template <typename T_result_t, typename c_layout_t, typename c_data_t>
80  static inline MFEM_ALWAYS_INLINE
81  void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
82  {
83  const int qpts = T_result_t::x_type::layout_type::dim_1;
84  const int ne = T_result_t::x_type::layout_type::dim_3;
85  const int vs = sizeof(T.x[0])/sizeof(T.x[0][0]);
86  for (int k = 0; k < ne; k++)
87  {
88  for (int i = 0; i < qpts; i++)
89  {
90  for (int s = 0; s < vs; s++)
91  {
92  c[l.ind(i,k)][s] = F.Eval1D(T.x(i,0,k)[s]);
93  }
94  }
95  }
96  }
97  };
98  template <bool dummy> struct Dim<2,dummy>
99  {
100  template <typename T_result_t, typename c_layout_t, typename c_data_t>
101  static inline MFEM_ALWAYS_INLINE
102  void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
103  {
104  const int qpts = T_result_t::x_type::layout_type::dim_1;
105  const int ne = T_result_t::x_type::layout_type::dim_3;
106  const int vs = sizeof(T.x[0])/sizeof(T.x[0][0]);
107  for (int k = 0; k < ne; k++)
108  {
109  for (int i = 0; i < qpts; i++)
110  {
111  for (int s = 0; s < vs; s++)
112  {
113  c[l.ind(i,k)][s] = F.Eval2D(T.x(i,0,k)[s], T.x(i,1,k)[s]);
114  }
115  }
116  }
117  }
118  };
119  template <bool dummy> struct Dim<3,dummy>
120  {
121  template <typename T_result_t, typename c_layout_t, typename c_data_t>
122  static inline MFEM_ALWAYS_INLINE
123  void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
124  {
125  const int qpts = T_result_t::x_type::layout_type::dim_1;
126  const int ne = T_result_t::x_type::layout_type::dim_3;
127  const int vs = sizeof(T.x[0])/sizeof(T.x[0][0]);
128  for (int k = 0; k < ne; k++)
129  {
130  for (int i = 0; i < qpts; i++)
131  {
132  for (int s = 0; s < vs; s++)
133  {
134  c[l.ind(i,k)][s] =
135  F.Eval3D(T.x(i,0,k)[s], T.x(i,1,k)[s], T.x(i,2,k)[s]);
136  }
137  }
138  }
139  }
140  };
141 
142 public:
143  /// Constructor for the case when Func has no data members.
145  /// Constructor for the case when Func has data members.
146  TFunctionCoefficient(Func &F_) : F(F_) { }
147  // Default copy constructor, Func has to have copy constructor.
148 
149  template <typename T_result_t, typename c_layout_t, typename c_data_t>
150  inline MFEM_ALWAYS_INLINE
151  void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
152  {
153  const int qpts = T_result_t::x_type::layout_type::dim_1;
154  const int sdim = T_result_t::x_type::layout_type::dim_2;
155  const int ne = T_result_t::x_type::layout_type::dim_3;
156  MFEM_STATIC_ASSERT(c_layout_t::rank == 2 && c_layout_t::dim_1 == qpts &&
157  c_layout_t::dim_2 == ne, "invalid c_layout_t");
158 
159  Dim<sdim,true>::Eval(F, T, l, c);
160  }
161 };
162 
163 
164 /// Piecewise constant coefficient class. The subdomains where the coefficient
165 /// is constant are given by the mesh attributes.
166 template <typename complex_t = double>
168 {
169 public:
170  static const bool uses_attributes = true;
171  typedef complex_t complex_type;
172 
173 protected:
174  Vector constants; // complex_t = double
175 
176 public:
177  /// Note: in the input array index i corresponds to mesh attribute i+1.
179  : constants(constants) { }
180  // default copy constructor
181 
182  template <typename T_result_t, typename c_layout_t, typename c_data_t>
183  inline MFEM_ALWAYS_INLINE
184  void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
185  {
186  const int ne = T_result_t::ne;
187  const int vs = sizeof(T.attrib[0])/sizeof(T.attrib[0][0]);
188  MFEM_STATIC_ASSERT(vs == sizeof(c[0])/sizeof(c[0][0]), "");
189  for (int i = 0; i < ne; i++)
190  {
191  typename c_data_t::data_type ci;
192  for (int s = 0; s < vs; s++)
193  {
194  ci[s] = constants(T.attrib[i][s]-1);
195  }
196  TAssign<AssignOp::Set>(l.ind2(i), c, ci);
197  }
198  }
199 };
200 
201 /// GridFunction coefficient class.
202 template <typename FieldEval>
204 {
205 public:
206  static const bool uses_element_idxs = true;
207 
208  typedef typename FieldEval::FESpace_type FESpace_type;
209  typedef typename FieldEval::ShapeEval_type ShapeEval_type;
210  typedef typename FieldEval::VecLayout_type VecLayout_type;
211  typedef typename FieldEval::complex_type complex_type;
212 
213 protected:
214  FieldEval fieldEval;
215 
216 public:
217  // This constructor uses a shallow copy of fE.fespace as part of fieldEval.
218  inline MFEM_ALWAYS_INLINE
219  TGridFunctionCoefficient(const FieldEval &fE,
220  const complex_type *data)
221  : fieldEval(fE, data, NULL)
222  { }
223 
224  // This constructor uses a shallow copy of tfes as part of fieldEval.
225  inline MFEM_ALWAYS_INLINE
227  const ShapeEval_type &shapeEval,
228  const VecLayout_type &vec_layout,
229  const complex_type *data)
230  : fieldEval(tfes, shapeEval, vec_layout, data, NULL)
231  { }
232 
233  // This constructor creates new FESpace_type as part of fieldEval.
234  inline MFEM_ALWAYS_INLINE
236  const complex_type *data)
237  : fieldEval(fes, data, NULL)
238  { }
239 
240  // This constructor creates new FESpace_type as part of fieldEval.
241  inline MFEM_ALWAYS_INLINE
243  : fieldEval(*func.FESpace(), func.GetData(), NULL)
244  { }
245 
246  // default copy constructor
247 
248  template <typename T_result_t, typename c_layout_t, typename c_data_t>
249  inline MFEM_ALWAYS_INLINE
250  void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
251  {
252  const int ne = T_result_t::ne;
253  const int vdim = FieldEval::vdim;
254  const int qpts = FieldEval::qpts;
255  MFEM_STATIC_ASSERT(c_layout_t::rank == 2, "tensor rank must be 2");
256  MFEM_STATIC_ASSERT(c_layout_t::dim_1 == qpts, "incompatible quadrature");
257  MFEM_STATIC_ASSERT(c_layout_t::dim_2 == ne, "");
258  MFEM_STATIC_ASSERT(vdim == 1, "vdim != 1 is not supported");
259 
260  fieldEval.GetValues(T.first_elem_idx, l.template split_2<1,ne>(), c);
261  }
262 };
263 
264 
265 /// Auxiliary class that is used to simplify the evaluation of a coefficient and
266 /// scaling it by the weights of a quadrature rule.
267 template <typename IR, typename coeff_t, typename impl_traits_t>
269 {
270  static const int qpts = IR::qpts;
271  static const int ne = impl_traits_t::batch_size;
272  typedef typename coeff_t::complex_type complex_type;
273  typedef typename impl_traits_t::vcomplex_t vcomplex_t;
274 
275  template <bool is_const, bool dummy> struct Aux;
276 
277  // constant coefficient
278  template <bool dummy> struct Aux<true,dummy>
279  {
280  typedef struct { } result_t;
282 
283  inline MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
284  {
285  c.Eval(true, cw.layout, cw);
286  int_rule.ApplyWeights(cw);
287  }
288 
289  template <typename T_result_t>
290  inline MFEM_ALWAYS_INLINE
291  void Eval(const T_result_t &F, result_t &res) { }
292 
293  inline MFEM_ALWAYS_INLINE
294  const complex_type &get(const result_t &res, int i, int k) const
295  {
296  return cw(i,0);
297  }
298  };
299  // non-constant coefficient
300  template <bool dummy> struct Aux<false,dummy>
301  {
303 #ifdef MFEM_TEMPLATE_INTRULE_COEFF_PRECOMP
305 #else
307 #endif
308  coeff_t c;
309 
310 #ifdef MFEM_TEMPLATE_INTRULE_COEFF_PRECOMP
311  inline MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
312  : c(c)
313  {
314  int_rule.template AssignWeights<AssignOp::Set>(w.layout, w);
315  }
316 #else
317  inline MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
318  : int_rule(int_rule), c(c) { }
319 #endif
320  template <typename T_result_t>
321  inline MFEM_ALWAYS_INLINE
322  void Eval(const T_result_t &F, result_t &res)
323  {
324  c.Eval(F, res.layout, res);
325 #ifdef MFEM_TEMPLATE_INTRULE_COEFF_PRECOMP
326  for (int i = 0; i < ne; i++)
327  {
328  TAssign<AssignOp::Mult>(res.layout.ind2(i), res,
329  w.layout.merge_12(), w);
330  }
331 #else
332  int_rule.template AssignWeights<AssignOp::Mult>(res.layout, res);
333 #endif
334  }
335 
336  inline MFEM_ALWAYS_INLINE
337  const vcomplex_t &get(const result_t &res, int i, int k) const
338  {
339  return res(i,k);
340  }
341  };
342 
344 };
345 
346 } // namespace mfem
347 
348 #endif // MFEM_TEMPLATE_COEFFICIENT
static const bool uses_Jacobians
TFunctionCoefficient(Func &F_)
Constructor for the case when Func has data members.
MFEM_ALWAYS_INLINE void Eval(const T_result_t &F, result_t &res)
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
Templated coefficient classes, cf. coefficient.?pp.
FieldEval::FESpace_type FESpace_type
TMatrix< qpts, ne, vcomplex_t > result_t
GridFunction coefficient class.
static const layout_type layout
Definition: ttensor.hpp:356
static const bool uses_coordinates
MFEM_ALWAYS_INLINE void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c) const
static MFEM_ALWAYS_INLINE void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
MFEM_ALWAYS_INLINE TGridFunctionCoefficient(const GridFunction &func)
FieldEval::complex_type complex_type
Function coefficient.
static const bool uses_element_idxs
TFunctionCoefficient()
Constructor for the case when Func has no data members.
FieldEval::ShapeEval_type ShapeEval_type
static const bool uses_element_idxs
static const bool uses_coordinates
static const bool is_const
MFEM_ALWAYS_INLINE void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
impl_traits_t::vcomplex_t vcomplex_t
TMatrix< qpts, 1, complex_type > cw
MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
static MFEM_ALWAYS_INLINE void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
MFEM_ALWAYS_INLINE TGridFunctionCoefficient(const FESpace_type &tfes, const ShapeEval_type &shapeEval, const VecLayout_type &vec_layout, const complex_type *data)
TPiecewiseConstCoefficient(const Vector &constants)
Note: in the input array index i corresponds to mesh attribute i+1.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
static const bool is_const
TConstantCoefficient(complex_t val)
TMatrix< qpts, 1, typename IR::real_type > w
FieldEval::VecLayout_type VecLayout_type
Aux< coeff_t::is_const, true > Type
MFEM_ALWAYS_INLINE void Eval(const T_result_t &F, result_t &res)
static MFEM_ALWAYS_INLINE void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
static OffsetStridedLayout1D< N1, S1 > ind2(int i2)
Definition: tlayout.hpp:115
Vector data type.
Definition: vector.hpp:58
MFEM_ALWAYS_INLINE void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
coeff_t::complex_type complex_type
RefCoord s[3]
MFEM_ALWAYS_INLINE TGridFunctionCoefficient(const FieldEval &fE, const complex_type *data)
static const int rank
static const bool uses_attributes
MFEM_ALWAYS_INLINE TGridFunctionCoefficient(const FiniteElementSpace &fes, const complex_type *data)
MFEM_ALWAYS_INLINE void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)