MFEM  v3.3
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tcoefficient.hpp
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 #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 // Function coefficient. The template class 'Func' has to implement at least one
60 // of the following methods, depending on the dimension that will be used:
61 // complex_t Eval1D(real_t);
62 // complex_t Eval2D(real_t,real_t);
63 // complex_t Eval3D(real_t,real_t,real_t);
64 // Use MFEM_FLOPS_ADD() to count flops inside Eval*D.
65 template <typename Func, typename complex_t = double>
67 {
68 public:
69  static const bool uses_coordinates = true;
70  typedef complex_t complex_type;
71 
72 protected:
73  Func F;
74 
75  template <int dim, bool dummy> struct Dim;
76  template <bool dummy> struct Dim<1,dummy>
77  {
78  template <typename T_result_t, typename c_layout_t, typename c_data_t>
79  static inline MFEM_ALWAYS_INLINE
80  void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
81  {
82  const int qpts = T_result_t::x_type::layout_type::dim_1;
83  const int ne = T_result_t::x_type::layout_type::dim_3;
84  for (int k = 0; k < ne; k++)
85  {
86  for (int i = 0; i < qpts; i++)
87  {
88  c[l.ind(i,k)] = F.Eval1D(T.x(i,0,k));
89  }
90  }
91  }
92  };
93  template <bool dummy> struct Dim<2,dummy>
94  {
95  template <typename T_result_t, typename c_layout_t, typename c_data_t>
96  static inline MFEM_ALWAYS_INLINE
97  void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
98  {
99  const int qpts = T_result_t::x_type::layout_type::dim_1;
100  const int ne = T_result_t::x_type::layout_type::dim_3;
101  for (int k = 0; k < ne; k++)
102  {
103  for (int i = 0; i < qpts; i++)
104  {
105  c[l.ind(i,k)] = F.Eval2D(T.x(i,0,k), T.x(i,1,k));
106  }
107  }
108  }
109  };
110  template <bool dummy> struct Dim<3,dummy>
111  {
112  template <typename T_result_t, typename c_layout_t, typename c_data_t>
113  static inline MFEM_ALWAYS_INLINE
114  void Eval(Func &F, const T_result_t &T, const c_layout_t &l, c_data_t &c)
115  {
116  const int qpts = T_result_t::x_type::layout_type::dim_1;
117  const int ne = T_result_t::x_type::layout_type::dim_3;
118  for (int k = 0; k < ne; k++)
119  {
120  for (int i = 0; i < qpts; i++)
121  {
122  c[l.ind(i,k)] = F.Eval3D(T.x(i,0,k), T.x(i,1,k), T.x(i,2,k));
123  }
124  }
125  }
126  };
127 
128 public:
129  // Constructor for the case when Func has no data members.
131  // Constructor for the case when Func has data members.
132  TFunctionCoefficient(Func &F_) : F(F_) { }
133  // Default copy constructor, Func has to have copy constructor.
134 
135  template <typename T_result_t, typename c_layout_t, typename c_data_t>
136  inline MFEM_ALWAYS_INLINE
137  void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
138  {
139  const int qpts = T_result_t::x_type::layout_type::dim_1;
140  const int sdim = T_result_t::x_type::layout_type::dim_2;
141  const int ne = T_result_t::x_type::layout_type::dim_3;
142  MFEM_STATIC_ASSERT(c_layout_t::rank == 2 && c_layout_t::dim_1 == qpts &&
143  c_layout_t::dim_2 == ne, "invalid c_layout_t");
144 
145  Dim<sdim,true>::Eval(F, T, l, c);
146  }
147 };
148 
149 
150 /// Piecewise constant coefficient class. The subdomains where the coefficient
151 /// is constant are given by the mesh attributes.
152 template <typename complex_t = double>
154 {
155 public:
156  static const bool uses_attributes = true;
157  typedef complex_t complex_type;
158 
159 protected:
160  Vector constants; // complex_t = double
161 
162 public:
163  /// Note: in the input array index i corresponds to mesh attribute i+1.
165  : constants(constants) { }
166  // default copy constructor
167 
168  template <typename T_result_t, typename c_layout_t, typename c_data_t>
169  inline MFEM_ALWAYS_INLINE
170  void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
171  {
172  const int ne = T_result_t::ne;
173  for (int i = 0; i < ne; i++)
174  {
175  TAssign<AssignOp::Set>(l.ind2(i), c, constants(T.attrib[i]-1));
176  }
177  }
178 };
179 
180 
181 template <typename FieldEval>
183 {
184 public:
185  static const bool uses_element_idxs = true;
186 
187  typedef typename FieldEval::FESpace_type FESpace_type;
188  typedef typename FieldEval::ShapeEval_type ShapeEval_type;
189  typedef typename FieldEval::VecLayout_type VecLayout_type;
190  typedef typename FieldEval::complex_type complex_type;
191 
192 protected:
193  FieldEval fieldEval;
194 
195 public:
196  // This constructor uses a shallow copy of fE.fespace as part of fieldEval.
197  inline MFEM_ALWAYS_INLINE
198  TGridFunctionCoefficient(const FieldEval &fE,
199  const complex_type *data)
200  : fieldEval(fE, data, NULL)
201  { }
202 
203  // This constructor uses a shallow copy of tfes as part of fieldEval.
204  inline MFEM_ALWAYS_INLINE
206  const ShapeEval_type &shapeEval,
207  const VecLayout_type &vec_layout,
208  const complex_type *data)
209  : fieldEval(tfes, shapeEval, vec_layout, data, NULL)
210  { }
211 
212  // This constructor creates new FESpace_type as part of fieldEval.
213  inline MFEM_ALWAYS_INLINE
215  const complex_type *data)
216  : fieldEval(fes, data, NULL)
217  { }
218 
219  // This constructor creates new FESpace_type as part of fieldEval.
220  inline MFEM_ALWAYS_INLINE
222  : fieldEval(*func.FESpace(), func.GetData(), NULL)
223  { }
224 
225  // default copy constructor
226 
227  template <typename T_result_t, typename c_layout_t, typename c_data_t>
228  inline MFEM_ALWAYS_INLINE
229  void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
230  {
231  const int ne = T_result_t::ne;
232  const int vdim = FieldEval::vdim;
233  const int qpts = FieldEval::qpts;
234  MFEM_STATIC_ASSERT(c_layout_t::rank == 2, "tensor rank must be 2");
235  MFEM_STATIC_ASSERT(c_layout_t::dim_1 == qpts, "incompatible quadrature");
236  MFEM_STATIC_ASSERT(c_layout_t::dim_2 == ne, "");
237  MFEM_STATIC_ASSERT(vdim == 1, "vdim != 1 is not supported");
238 
239  fieldEval.GetValues(T.first_elem_idx, l.template split_2<1,ne>(), c);
240  }
241 };
242 
243 
244 /// Auxiliary class that is used to simplify the evaluation of a coefficient and
245 /// scaling it by the weights of a quadrature rule.
246 template <typename IR, typename coeff_t, int NE>
248 {
249  static const int qpts = IR::qpts;
250  static const int ne = NE;
252 
253  template <bool is_const, bool dummy> struct Aux;
254 
255  // constant coefficient
256  template <bool dummy> struct Aux<true,dummy>
257  {
258  typedef struct { } result_t;
260 
261  inline MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
262  {
263  c.Eval(true, cw.layout, cw);
264  int_rule.ApplyWeights(cw);
265  }
266 
267  template <typename T_result_t>
268  inline MFEM_ALWAYS_INLINE
269  void Eval(const T_result_t &F, result_t &res) { }
270 
271  inline MFEM_ALWAYS_INLINE
272  const complex_type &get(const result_t &res, int i, int k) const
273  {
274  return cw(i,0);
275  }
276  };
277  // non-constant coefficient
278  template <bool dummy> struct Aux<false,dummy>
279  {
281 #ifdef MFEM_TEMPLATE_INTRULE_COEFF_PRECOMP
283 #else
285 #endif
287 
288 #ifdef MFEM_TEMPLATE_INTRULE_COEFF_PRECOMP
289  inline MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
290  : c(c)
291  {
292  int_rule.template AssignWeights<AssignOp::Set>(w.layout, w);
293  }
294 #else
295  inline MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
296  : int_rule(int_rule), c(c) { }
297 #endif
298  template <typename T_result_t>
299  inline MFEM_ALWAYS_INLINE
300  void Eval(const T_result_t &F, result_t &res)
301  {
302  c.Eval(F, res.layout, res);
303 #ifdef MFEM_TEMPLATE_INTRULE_COEFF_PRECOMP
304  for (int i = 0; i < ne; i++)
305  {
306  TAssign<AssignOp::Mult>(res.layout.ind2(i), res,
307  w.layout.merge_12(), w);
308  }
309 #else
310  int_rule.template AssignWeights<AssignOp::Mult>(res.layout, res);
311 #endif
312  }
313 
314  inline MFEM_ALWAYS_INLINE
315  const complex_type &get(const result_t &res, int i, int k) const
316  {
317  return res(i,k);
318  }
319  };
320 
322 };
323 
324 } // namespace mfem
325 
326 #endif // MFEM_TEMPLATE_COEFFICIENT
TMatrix< qpts, ne, complex_type > result_t
static const bool uses_Jacobians
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
FieldEval::FESpace_type FESpace_type
static const layout_type layout
Definition: ttensor.hpp:311
static const bool uses_coordinates
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 GridFunction &func)
FieldEval::complex_type complex_type
static const bool uses_element_idxs
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)
coeff_t::complex_type complex_type
TConstantCoefficient coeff_t
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)
MFEM_ALWAYS_INLINE Aux(const IR &int_rule, const coeff_t &c)
Aux< coeff_t::is_const, true > Type
TPiecewiseConstCoefficient(const Vector &constants)
Note: in the input array index i corresponds to mesh attribute i+1.
MFEM_ALWAYS_INLINE void Eval(const T_result_t &F, result_t &res)
static StridedLayout1D< N1 *N2,(S1< S2)?S1:S2 > merge_12()
Definition: tlayout.hpp:149
static const bool is_const
TConstantCoefficient(complex_t val)
MFEM_ALWAYS_INLINE void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c) const
FieldEval::VecLayout_type VecLayout_type
MFEM_ALWAYS_INLINE void Eval(const T_result_t &F, result_t &res)
TMatrix< qpts, 1, complex_type > cw
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:114
TMatrix< qpts, 1, typename IR::real_type > w
Vector data type.
Definition: vector.hpp:36
MFEM_ALWAYS_INLINE void Eval(const T_result_t &T, const c_layout_t &l, c_data_t &c)
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)