MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tbilinearform.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_BILINEAR_FORM
13 #define MFEM_TEMPLATE_BILINEAR_FORM
14 
15 #include "../config/tconfig.hpp"
16 #include "../linalg/ttensor.hpp"
17 #include "bilinearform.hpp"
18 #include "tevaluator.hpp"
19 #include "teltrans.hpp"
20 #include "tcoefficient.hpp"
21 #include "fespace.hpp"
22 
23 namespace mfem
24 {
25 
26 // Templated bilinear form class, cf. bilinearform.?pp
27 
28 // complex_t - sol dof data type
29 // real_t - mesh nodes, sol basis, mesh basis data type
30 template <typename meshType, typename solFESpace,
31  typename IR, typename IntegratorType,
32  typename solVecLayout_t = ScalarLayout,
33  typename complex_t = double, typename real_t = double>
34 class TBilinearForm : public Operator
35 {
36 protected:
37  typedef complex_t complex_type;
38  typedef real_t real_type;
39 
40  typedef typename meshType::FE_type meshFE_type;
42  typedef typename solFESpace::FE_type solFE_type;
44  typedef solVecLayout_t solVecLayout_type;
45 
46  static const int dim = meshType::dim;
47  static const int sdim = meshType::space_dim;
48  static const int dofs = solFE_type::dofs;
49  static const int vdim = solVecLayout_t::vec_dim;
50  static const int qpts = IR::qpts;
51 
52  typedef IntegratorType integ_t;
53  typedef typename integ_t::coefficient_type coeff_t;
54  typedef typename integ_t::template kernel<sdim,dim,complex_t>::type kernel_t;
55  typedef typename kernel_t::template p_asm_data<qpts>::type p_assembled_t;
56  typedef typename kernel_t::template f_asm_data<qpts>::type f_assembled_t;
57 
59  template <int NE> struct T_result
60  {
61  static const int EvalOps =
62  Trans_t::template Get<coeff_t,kernel_t>::EvalOps;
63  typedef typename Trans_t::template Result<EvalOps,NE> Type;
64  };
65 
66  typedef FieldEvaluator<solFESpace,solVecLayout_t,IR,
67  complex_t,real_t> solFieldEval;
68  template <int BE> struct S_spec
69  {
70  typedef typename solFieldEval::template Spec<kernel_t,BE> Spec;
71  typedef typename Spec::DataType DataType;
72  typedef typename Spec::ElementMatrix ElementMatrix;
73  };
74 
75  // Data members
76 
77  meshType mesh;
79 
82  mutable solFESpace solFES;
83  solVecLayout_t solVecLayout;
84 
86 
88 
90 
92 
93 public:
94  TBilinearForm(const IntegratorType &integ, const FiniteElementSpace &sol_fes)
95  : Operator(sol_fes.GetNDofs()*vdim),
96  mesh(*sol_fes.GetMesh()),
97  meshEval(mesh.fe),
98  sol_fe(*sol_fes.FEColl()),
99  solEval(sol_fe),
100  solFES(sol_fe, sol_fes),
101  solVecLayout(sol_fes),
102  int_rule(),
103  coeff(integ.coeff),
104  assembled_data(NULL),
105  in_fes(sol_fes)
106  { }
107 
108  virtual ~TBilinearForm()
109  {
110  delete [] assembled_data;
111  }
112 
113  /// Get the input finite element space prolongation matrix
114  virtual const Operator *GetProlongation() const
115  { return ((FiniteElementSpace &)in_fes).GetProlongationMatrix(); }
116  /// Get the input finite element space restriction matrix
117  virtual const Operator *GetRestriction() const
118  { return ((FiniteElementSpace &)in_fes).GetRestrictionMatrix(); }
119 
120  virtual void Mult(const Vector &x, Vector &y) const
121  {
122  if (assembled_data)
123  {
124  const int num_elem = 1;
125  MultAssembled<num_elem>(x, y);
126  }
127  else
128  {
129  MultUnassembled(x, y);
130  }
131  }
132 
133  // complex_t = double
134  void MultUnassembled(const Vector &x, Vector &y) const
135  {
136  y = 0.0;
137 
138  const int BE = 1; // batch-size of elements
139  typedef typename kernel_t::template
140  CoefficientEval<IR,coeff_t,BE>::Type coeff_eval_t;
141 
142  // For better performance, create stack copies of solFES, and solEval
143  // inside 'solFEval'. The element-transformation 'T' also copies the
144  // meshFES, meshEval, etc internally.
145  // Is performance actually better with this implementation?
146  Trans_t T(mesh, meshEval);
148  x.GetData(), y.GetData());
149  coeff_eval_t wQ(int_rule, coeff);
150 
151  const int NE = mesh.GetNE();
152  for (int el = 0; el < NE; el++)
153  {
154 #if 0
155  typename S_spec<BE>::DataType R;
156  solFEval.Eval(el, R);
157 
158  typename T_result<BE>::Type F;
159  T.Eval(el, F);
160 #else
161  typename T_result<BE>::Type F;
162  T.Eval(el, F);
163 
164  typename S_spec<BE>::DataType R;
165  solFEval.Eval(el, R);
166 #endif
167 
168  typename coeff_eval_t::result_t res;
169  wQ.Eval(F, res);
170 
171  kernel_t::Action(0, F, wQ, res, R);
172 
173  solFEval.template Assemble<true>(R);
174  }
175  }
176 
177  // Partial assembly of quadrature point data
178  void Assemble()
179  {
180  const int BE = 1; // batch-size of elements
181  typedef typename kernel_t::template
182  CoefficientEval<IR,coeff_t,BE>::Type coeff_eval_t;
183 
184  Trans_t T(mesh, meshEval);
185  coeff_eval_t wQ(int_rule, coeff);
186 
187  const int NE = mesh.GetNE();
188  if (!assembled_data)
189  {
190  assembled_data = new p_assembled_t[NE];
191  }
192  for (int el = 0; el < NE; el++) // BE == 1
193  {
194  typename T_result<BE>::Type F;
195  T.Eval(el, F);
196 
197  typename coeff_eval_t::result_t res;
198  wQ.Eval(F, res);
199 
200  for (int k = 0; k < BE; k++)
201  {
202  kernel_t::Assemble(k, F, wQ, res, assembled_data[el+k]);
203  }
204  }
205  }
206 
207  template <int num_elem>
208  inline MFEM_ALWAYS_INLINE
209  void ElementAddMultAssembled(int el, solFieldEval &solFEval) const
210  {
211  typename S_spec<num_elem>::DataType R;
212  solFEval.Eval(el, R);
213 
214  for (int k = 0; k < num_elem; k++)
215  {
216  kernel_t::MultAssembled(k, assembled_data[el+k], R);
217  }
218 
219  solFEval.template Assemble<true>(R);
220  }
221 
222  // complex_t = double
223  template <int num_elem>
224  void MultAssembled(const Vector &x, Vector &y) const
225  {
226  y = 0.0;
227 
229  x.GetData(), y.GetData());
230 
231  const int NE = mesh.GetNE();
232  const int bNE = NE-NE%num_elem;
233  for (int el = 0; el < bNE; el += num_elem)
234  {
235  ElementAddMultAssembled<num_elem>(el, solFEval);
236  }
237  for (int el = bNE; el < NE; el++)
238  {
239  ElementAddMultAssembled<1>(el, solFEval);
240  }
241  }
242 
243 #ifdef MFEM_TEMPLATE_ENABLE_SERIALIZE
244  // complex_t = double
246  {
247  y = 0.0;
248 
250  solFESpace solFES(this->solFES);
251 
253 
254  const int NE = mesh.GetNE();
255  for (int el = 0; el < NE; el++)
256  {
257  solFES.SetElement(el);
258 
259  solFES.VectorExtract(solVecLayout, x, xy_dof.layout, xy_dof);
260  solFES.VectorAssemble(xy_dof.layout, xy_dof, solVecLayout, y);
261  }
262  }
263 
264  // real_t = double
265  void SerializeNodes(Vector &sNodes) const
266  {
267  typedef typename meshType::FESpace_type meshFESpace;
268  meshFESpace meshFES(mesh.t_fes);
270 
271  const int NE = mesh.GetNE();
272  sNodes.SetSize(lnodes_t::size*NE);
273  real_t *lNodes = sNodes.GetData();
274  for (int el = 0; el < NE; el++)
275  {
276  meshFES.SetElement(el);
277  meshFES.VectorExtract(mesh.node_layout, mesh.Nodes,
278  lnodes_t::layout, lNodes);
279  lNodes += lnodes_t::size;
280  }
281  }
282 
283  // partial assembly from "serialized" nodes
284  // real_t = double
286  {
287  const int BE = 1; // batch-size of elements
288  typedef typename kernel_t::template
289  CoefficientEval<IR,coeff_t,BE>::Type coeff_eval_t;
290 
291  Trans_t T(this->mesh, this->meshEval);
292  coeff_eval_t wQ(int_rule, coeff);
293 
294  const int NE = mesh.GetNE();
295  if (!assembled_data)
296  {
297  assembled_data = new p_assembled_t[NE];
298  }
299  for (int el = 0; el < NE; el++)
300  {
301  typename T_result<BE>::Type F;
302  T.EvalSerialized(el, sNodes.GetData(), F);
303 
304  typename coeff_eval_t::result_t res;
305  wQ.Eval(F, res);
306 
307  kernel_t::Assemble(0, F, wQ, res, assembled_data[el]);
308  }
309  }
310 
311  // complex_t = double
312  void Serialize(const Vector &x, Vector &sx) const
313  {
314  solVecLayout_t solVecLayout(this->solVecLayout);
315  typedef TTensor3<dofs,vdim,1,complex_t> vdof_data_t;
316  solFESpace solFES(this->solFES);
317 
318  const int NE = mesh.GetNE();
319  sx.SetSize(vdim*dofs*NE);
320  complex_t *loc_sx = sx.GetData();
321  for (int el = 0; el < NE; el++)
322  {
323  solFES.SetElement(el);
324  solFES.VectorExtract(solVecLayout, x, vdof_data_t::layout, loc_sx);
325  loc_sx += vdim*dofs;
326  }
327  }
328 
329  // serialized vector sx --> serialized vector 'sy'
330  // complex_t = double
331  void MultAssembledSerialized(const Vector &sx, Vector &sy) const
332  {
333  solFieldEval solFEval(solFES, solEval, solVecLayout, NULL, NULL);
334 
335  const int NE = mesh.GetNE();
336  const complex_t *loc_sx = sx.GetData();
337  complex_t *loc_sy = sy.GetData();
338  for (int el = 0; el < NE; el++)
339  {
340  typename S_spec<1>::DataType R;
341  solFEval.EvalSerialized(loc_sx, R);
342 
343  kernel_t::MultAssembled(0, assembled_data[el], R);
344 
345  solFEval.template AssembleSerialized<false>(R, loc_sy);
346 
347  loc_sx += vdim*dofs;
348  loc_sy += vdim*dofs;
349  }
350  }
351 #endif // MFEM_TEMPLATE_ENABLE_SERIALIZE
352 
353  // Assemble the operator in a SparseMatrix.
354  // complex_t = double
356  {
357  const int BE = 1; // batch-size of elements
358  typedef typename kernel_t::template
359  CoefficientEval<IR,coeff_t,BE>::Type coeff_eval_t;
360 
361  Trans_t T(mesh, meshEval);
362  solFESpace solFES(this->solFES);
364  solVecLayout_t solVecLayout(this->solVecLayout);
365  coeff_eval_t wQ(int_rule, coeff);
366 
367  const int NE = mesh.GetNE();
368  for (int el = 0; el < NE; el++)
369  {
370  f_assembled_t asm_qpt_data;
371  {
372  typename T_result<BE>::Type F;
373  T.Eval(el, F);
374 
375  typename coeff_eval_t::result_t res;
376  wQ.Eval(F, res);
377 
378  kernel_t::Assemble(0, F, wQ, res, asm_qpt_data);
379  }
380 
381  // For now, when vdim > 1, assume block-diagonal matrix with the same
382  // diagonal block for all components.
383  TMatrix<dofs,dofs> M_loc;
385  asm_qpt_data.layout, asm_qpt_data, M_loc.layout, M_loc, solEval);
386 
387  solFES.SetElement(el);
388  for (int bi = 0; bi < vdim; bi++)
389  {
390  solFES.AssembleBlock(bi, bi, solVecLayout, M_loc, M);
391  }
392  }
393  }
394 
395  // Assemble element matrices and store them as a DenseTensor object.
396  // complex_t = double
398  {
399  const int BE = 1; // batch-size of elements
400  typedef typename kernel_t::template
401  CoefficientEval<IR,coeff_t,BE>::Type coeff_eval_t;
402 
403  Trans_t T(mesh, meshEval);
405  coeff_eval_t wQ(int_rule, coeff);
406 
407  const int NE = mesh.GetNE();
408  for (int el = 0; el < NE; el++)
409  {
410  f_assembled_t asm_qpt_data;
411  {
412  typename T_result<BE>::Type F;
413  T.Eval(el, F);
414 
415  typename coeff_eval_t::result_t res;
416  wQ.Eval(F, res);
417 
418  kernel_t::Assemble(0, F, wQ, res, asm_qpt_data);
419  }
420 
421  // For now, when vdim > 1, assume block-diagonal matrix with the same
422  // diagonal block for all components.
423  // M is assumed to be (dof x dof x NE).
424  TMatrix<dofs,dofs> M_loc;
426  asm_qpt_data.layout, asm_qpt_data, M_loc.layout, M_loc, solEval);
427 
428  complex_t *M_data = M.GetData(el);
429  M_loc.template AssignTo<AssignOp::Set>(M_data);
430  }
431  }
432 
433  // Assemble element matrices and add them to the bilinear form
434  // complex_t = double
436  {
437  const int BE = 1; // batch-size of elements
438  typedef typename kernel_t::template
439  CoefficientEval<IR,coeff_t,BE>::Type coeff_eval_t;
440 
441  Trans_t T(mesh, meshEval);
443  coeff_eval_t wQ(int_rule, coeff);
444 
445  Array<int> vdofs;
446  const Array<int> *dof_map = sol_fe.GetDofMap();
447  const int *dof_map_ = dof_map->GetData();
448  DenseMatrix M_loc_perm(dofs*vdim,dofs*vdim); // initialized with zeros
449 
450  const int NE = mesh.GetNE();
451  for (int el = 0; el < NE; el++)
452  {
453  f_assembled_t asm_qpt_data;
454  {
455  typename T_result<BE>::Type F;
456  T.Eval(el, F);
457 
458  typename coeff_eval_t::result_t res;
459  wQ.Eval(F, res);
460 
461  kernel_t::Assemble(0, F, wQ, res, asm_qpt_data);
462  }
463 
464  // For now, when vdim > 1, assume block-diagonal matrix with the same
465  // diagonal block for all components.
466  TMatrix<dofs,dofs> M_loc;
468  asm_qpt_data.layout, asm_qpt_data, M_loc.layout, M_loc, solEval);
469 
470  if (dof_map) // switch from tensor-product ordering
471  {
472  for (int i = 0; i < dofs; i++)
473  {
474  for (int j = 0; j < dofs; j++)
475  {
476  M_loc_perm(dof_map_[i],dof_map_[j]) = M_loc(i,j);
477  }
478  }
479  for (int bi = 1; bi < vdim; bi++)
480  {
481  M_loc_perm.CopyMN(M_loc_perm, dofs, dofs, 0, 0,
482  bi*dofs, bi*dofs);
483  }
484  a.AssembleElementMatrix(el, M_loc_perm, vdofs);
485  }
486  else
487  {
488  DenseMatrix DM(M_loc.data, dofs, dofs);
489  if (vdim == 1)
490  {
491  a.AssembleElementMatrix(el, DM, vdofs);
492  }
493  else
494  {
495  for (int bi = 0; bi < vdim; bi++)
496  {
497  M_loc_perm.CopyMN(DM, dofs, dofs, 0, 0, bi*dofs, bi*dofs);
498  }
499  a.AssembleElementMatrix(el, M_loc_perm, vdofs);
500  }
501  }
502  }
503  }
504 
505  // Multiplication using assembled element matrices stored as a DenseTensor.
506  // complex_t = double
507  void AddMult(DenseTensor &M, const Vector &x, Vector &y) const
508  {
509  // For now, when vdim > 1, assume block-diagonal matrix with the same
510  // diagonal block for all components.
511  // M is assumed to be (dof x dof x NE).
512  solVecLayout_t solVecLayout(this->solVecLayout);
513  const int NE = mesh.GetNE();
514  for (int el = 0; el < NE; el++)
515  {
516  TTensor3<dofs,vdim,1,complex_t> x_dof, y_dof;
517 
518  solFES.SetElement(el);
519  solFES.VectorExtract(solVecLayout, x, x_dof.layout, x_dof);
520  Mult_AB<false>(TMatrix<dofs,dofs>::layout,
521  M(el).Data(),
522  x_dof.layout.merge_23(), x_dof,
523  y_dof.layout.merge_23(), y_dof);
524  solFES.VectorAssemble(y_dof.layout, y_dof, solVecLayout, y);
525  }
526  }
527 };
528 
529 } // namespace mfem
530 
531 #endif // MFEM_TEMPLATE_BILINEAR_FORM
meshType::FE_type meshFE_type
FieldEvaluator< solFESpace, solVecLayout_t, IR, complex_t, real_t > solFieldEval
MFEM_ALWAYS_INLINE void EvalSerialized(int el, const real_t *nodeData, Result< EvalOps, NE > &F)
Definition: teltrans.hpp:119
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
static const layout_type layout
Definition: ttensor.hpp:311
const FiniteElementSpace & in_fes
solVecLayout_t solVecLayout
MFEM_ALWAYS_INLINE void Eval(DataType &F)
T * GetData()
Returns the data.
Definition: array.hpp:92
integ_t::coefficient_type coeff_t
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void AssembleBilinearForm(BilinearForm &a) const
kernel_t::template p_asm_data< qpts >::type p_assembled_t
void AddMult(DenseTensor &M, const Vector &x, Vector &y) const
virtual const Operator * GetProlongation() const
Get the input finite element space prolongation matrix.
static const layout_type layout
Definition: ttensor.hpp:347
double * GetData() const
Definition: vector.hpp:121
void AssembleFromSerializedNodes(const Vector &sNodes)
static const int dim
void MultUnassembled(const Vector &x, Vector &y) const
int dim
Definition: ex3.cpp:47
TBilinearForm(const IntegratorType &integ, const FiniteElementSpace &sol_fes)
Data type sparse matrix.
Definition: sparsemat.hpp:38
MFEM_ALWAYS_INLINE void EvalSerialized(const complex_t *loc_dofs, DataType &F)
solFESpace::FE_type solFE_type
double * GetData(int k)
Definition: densemat.hpp:692
void Serialize(const Vector &x, Vector &sx) const
virtual const Operator * GetRestriction() const
Get the input finite element space restriction matrix.
Trans_t::template Result< EvalOps, NE > Type
meshShapeEval meshEval
solVecLayout_t solVecLayout_type
ShapeEvaluator< solFE_type, IR, real_t > solShapeEval
void AssembleElementMatrix(int i, const DenseMatrix &elmat, Array< int > &vdofs, int skip_zeros=1)
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void TestElementwiseExtractAssemble(const Vector &x, Vector &y) const
static const int dofs
MFEM_ALWAYS_INLINE void ElementAddMultAssembled(int el, solFieldEval &solFEval) const
static const int sdim
MFEM_ALWAYS_INLINE void Eval(int el, Result< EvalOps, NE > &F)
Definition: teltrans.hpp:111
IntegratorType integ_t
TElementTransformation< meshType, IR, real_t > Trans_t
void MultAssembledSerialized(const Vector &sx, Vector &sy) const
solFieldEval::template Spec< kernel_t, BE > Spec
void AssembleMatrix(DenseTensor &M) const
integ_t::template kernel< sdim, dim, complex_t >::type kernel_t
static const int qpts
Vector data type.
Definition: vector.hpp:41
ShapeEvaluator< meshFE_type, IR, real_t > meshShapeEval
void SerializeNodes(Vector &sNodes) const
static const int vdim
p_assembled_t * assembled_data
void CopyMN(const DenseMatrix &A, int m, int n, int Aro, int Aco)
Copy the m x n submatrix of A at row/col offsets Aro/Aco to *this.
Definition: densemat.cpp:2567
double * Data()
Definition: densemat.hpp:694
void MultAssembled(const Vector &x, Vector &y) const
data_t data[aligned_size >0?aligned_size:1]
Definition: ttensor.hpp:249
solShapeEval solEval
kernel_t::template f_asm_data< qpts >::type f_assembled_t
Abstract operator.
Definition: operator.hpp:21
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:634
void AssembleMatrix(SparseMatrix &M) const
Spec::ElementMatrix ElementMatrix