MFEM  v3.3
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
tfespace.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_FESPACE
13 #define MFEM_TEMPLATE_FESPACE
14 
15 #include "../config/tconfig.hpp"
16 #include "../general/tassign.hpp"
17 #include "../linalg/ttensor.hpp"
18 #include "tfe.hpp" // for TFiniteElementSpace_simple
19 #include "fespace.hpp"
20 
21 namespace mfem
22 {
23 
24 // Templated finite element space classes, cf. fespace.?pp and fe_coll.?pp
25 
26 // Index types
27 
28 // IndexType must define:
29 // - constructor IndexType(const FE &fe, const FiniteElementSpace &fes)
30 // - copy constructor
31 // - void SetElement(int elem_idx)
32 // - int map(int loc_dof_idx, int elem_offset) const --> glob_dof_idx, for
33 // single component; elem_offset is relative to the currently set element.
34 
35 // Index type based on an array listing the dofs for each element where all
36 // elements are assumed to have the same number of dofs. Such an array is
37 // constructed from the J array of an element-to-dof Table with optional local
38 // renumbering to ensure tensor-product local dof ordering when needed.
39 template <typename FE>
41 {
42 protected:
43  const int *el_dof_list, *loc_dof_list;
44  bool own_list;
45 
46 public:
47  typedef FE FE_type;
48 
49  ElementDofIndexer(const FE &fe, const FiniteElementSpace &fes)
50  {
51  const Array<int> *loc_dof_map = fe.GetDofMap();
52  // fes.BuildElementToDofTable();
53  const Table &el_dof = fes.GetElementToDofTable();
54  MFEM_ASSERT(el_dof.Size_of_connections() == el_dof.Size() * FE::dofs,
55  "the element-to-dof Table is not compatible with this FE!");
56  int num_dofs = el_dof.Size() * FE::dofs;
57  if (!loc_dof_map)
58  {
59  // no local dof reordering
60  el_dof_list = el_dof.GetJ();
61  own_list = false;
62  }
63  else
64  {
65  // reorder the local dofs according to loc_dof_map
66  int *el_dof_list_ = new int[num_dofs];
67  const int *loc_dof_map_ = loc_dof_map->GetData();
68  for (int i = 0; i < el_dof.Size(); i++)
69  {
70  MFEM_ASSERT(el_dof.RowSize(i) == FE::dofs,
71  "incompatible element-to-dof Table!");
72  for (int j = 0; j < FE::dofs; j++)
73  {
74  el_dof_list_[j+FE::dofs*i] =
75  el_dof.GetJ()[loc_dof_map_[j]+FE::dofs*i];
76  }
77  }
78  el_dof_list = el_dof_list_;
79  own_list = true;
80  }
81  loc_dof_list = el_dof_list; // point to element 0
82  }
83 
84  // Shallow copy constructor
85  inline MFEM_ALWAYS_INLINE
87  : el_dof_list(orig.el_dof_list),
89  own_list(false)
90  { }
91 
92  inline MFEM_ALWAYS_INLINE
93  ~ElementDofIndexer() { if (own_list) { delete [] el_dof_list; } }
94 
95  inline MFEM_ALWAYS_INLINE
96  void SetElement(int elem_idx)
97  {
98  loc_dof_list = el_dof_list + elem_idx * FE::dofs;
99  }
100 
101  inline MFEM_ALWAYS_INLINE
102  int map(int loc_dof_idx, int elem_offset) const
103  {
104  return loc_dof_list[loc_dof_idx + elem_offset * FE::dofs];
105  }
106 };
107 
108 
109 // Simple template Finite Element Space, built using an IndexType. For a
110 // description of the requirements on IndexType, see above.
111 template <typename FE, typename IndexType>
113 {
114 public:
115  typedef FE FE_type;
116  typedef IndexType index_type;
117 
118 protected:
120 
121 public:
123  : ind(fe, fes) { }
124 
125  // default copy constructor
126 
127  void SetElement(int el) { ind.SetElement(el); }
128 
129  // Multi-element Extract:
130  // Extract dofs for multiple elements starting with the current element.
131  // The number of elements to extract is given by the second dimension of
132  // dof_layout_t: dof_layout is (DOFS x NumElems).
133  template <AssignOp::Type Op, typename glob_dof_data_t,
134  typename dof_layout_t, typename dof_data_t>
135  inline MFEM_ALWAYS_INLINE
136  void Extract(const glob_dof_data_t &glob_dof_data,
137  const dof_layout_t &dof_layout,
138  dof_data_t &dof_data) const
139  {
140  const int NE = dof_layout_t::dim_2;
141  MFEM_STATIC_ASSERT(FE::dofs == dof_layout_t::dim_1,
142  "invalid number of dofs");
143  for (int j = 0; j < NE; j++)
144  {
145  for (int i = 0; i < FE::dofs; i++)
146  {
147  Assign<Op>(dof_data[dof_layout.ind(i,j)],
148  glob_dof_data[ind.map(i,j)]);
149  }
150  }
151  }
152 
153  template <typename glob_dof_data_t,
154  typename dof_layout_t, typename dof_data_t>
155  inline MFEM_ALWAYS_INLINE
156  void Extract(const glob_dof_data_t &glob_dof_data,
157  const dof_layout_t &dof_layout,
158  dof_data_t &dof_data) const
159  {
160  Extract<AssignOp::Set>(glob_dof_data, dof_layout, dof_data);
161  }
162 
163  // Multi-element assemble.
164  template <AssignOp::Type Op,
165  typename dof_layout_t, typename dof_data_t,
166  typename glob_dof_data_t>
167  inline MFEM_ALWAYS_INLINE
168  void Assemble(const dof_layout_t &dof_layout,
169  const dof_data_t &dof_data,
170  glob_dof_data_t &glob_dof_data) const
171  {
172  const int NE = dof_layout_t::dim_2;
173  MFEM_STATIC_ASSERT(FE::dofs == dof_layout_t::dim_1,
174  "invalid number of dofs");
175  for (int j = 0; j < NE; j++)
176  {
177  for (int i = 0; i < FE::dofs; i++)
178  {
179  Assign<Op>(glob_dof_data[ind.map(i,j)],
180  dof_data[dof_layout.ind(i,j)]);
181  }
182  }
183  }
184 
185  template <typename dof_layout_t, typename dof_data_t,
186  typename glob_dof_data_t>
187  inline MFEM_ALWAYS_INLINE
188  void Assemble(const dof_layout_t &dof_layout,
189  const dof_data_t &dof_data,
190  glob_dof_data_t &glob_dof_data) const
191  {
192  Assemble<AssignOp::Add>(dof_layout, dof_data, glob_dof_data);
193  }
194 
195  // Multi-element VectorExtract: vdof_layout is (DOFS x NumComp x NumElems).
196  template <AssignOp::Type Op,
197  typename vec_layout_t, typename glob_vdof_data_t,
198  typename vdof_layout_t, typename vdof_data_t>
199  inline MFEM_ALWAYS_INLINE
200  void VectorExtract(const vec_layout_t &vl,
201  const glob_vdof_data_t &glob_vdof_data,
202  const vdof_layout_t &vdof_layout,
203  vdof_data_t &vdof_data) const
204  {
205  const int NC = vdof_layout_t::dim_2;
206  const int NE = vdof_layout_t::dim_3;
207  MFEM_STATIC_ASSERT(FE::dofs == vdof_layout_t::dim_1,
208  "invalid number of dofs");
209  MFEM_ASSERT(NC == vl.NumComponents(), "invalid number of components");
210  for (int k = 0; k < NC; k++)
211  {
212  for (int j = 0; j < NE; j++)
213  {
214  for (int i = 0; i < FE::dofs; i++)
215  {
216  Assign<Op>(vdof_data[vdof_layout.ind(i,k,j)],
217  glob_vdof_data[vl.ind(ind.map(i,j), k)]);
218  }
219  }
220  }
221  }
222 
223  template <typename vec_layout_t, typename glob_vdof_data_t,
224  typename vdof_layout_t, typename vdof_data_t>
225  inline MFEM_ALWAYS_INLINE
226  void VectorExtract(const vec_layout_t &vl,
227  const glob_vdof_data_t &glob_vdof_data,
228  const vdof_layout_t &vdof_layout,
229  vdof_data_t &vdof_data) const
230  {
231  VectorExtract<AssignOp::Set>(vl, glob_vdof_data, vdof_layout, vdof_data);
232  }
233 
234  // Multi-element VectorAssemble: vdof_layout is (DOFS x NumComp x NumElems).
235  template <AssignOp::Type Op,
236  typename vdof_layout_t, typename vdof_data_t,
237  typename vec_layout_t, typename glob_vdof_data_t>
238  inline MFEM_ALWAYS_INLINE
239  void VectorAssemble(const vdof_layout_t &vdof_layout,
240  const vdof_data_t &vdof_data,
241  const vec_layout_t &vl,
242  glob_vdof_data_t &glob_vdof_data) const
243  {
244  const int NC = vdof_layout_t::dim_2;
245  const int NE = vdof_layout_t::dim_3;
246  MFEM_STATIC_ASSERT(FE::dofs == vdof_layout_t::dim_1,
247  "invalid number of dofs");
248  MFEM_ASSERT(NC == vl.NumComponents(), "invalid number of components");
249  for (int k = 0; k < NC; k++)
250  {
251  for (int j = 0; j < NE; j++)
252  {
253  for (int i = 0; i < FE::dofs; i++)
254  {
255  Assign<Op>(glob_vdof_data[vl.ind(ind.map(i,j), k)],
256  vdof_data[vdof_layout.ind(i,k,j)]);
257  }
258  }
259  }
260  }
261 
262  template <typename vdof_layout_t, typename vdof_data_t,
263  typename vec_layout_t, typename glob_vdof_data_t>
264  inline MFEM_ALWAYS_INLINE
265  void VectorAssemble(const vdof_layout_t &vdof_layout,
266  const vdof_data_t &vdof_data,
267  const vec_layout_t &vl,
268  glob_vdof_data_t &glob_vdof_data) const
269  {
270  VectorAssemble<AssignOp::Add>(vdof_layout, vdof_data, vl, glob_vdof_data);
271  }
272 
273  // Extract a static number of consecutive components; vdof_layout is
274  // (dofs x NC x NE), where NC is the number of components to extract. It is
275  // assumed that: first_comp + NC <= vl.NumComponents().
276  template <typename vdof_layout_t, typename vdof_data_t,
277  typename vec_layout_t, typename glob_vdof_data_t>
278  inline MFEM_ALWAYS_INLINE
279  void ExtractComponents(int first_comp,
280  const vec_layout_t &vl,
281  const glob_vdof_data_t &glob_vdof_data,
282  const vdof_layout_t &vdof_layout,
283  vdof_data_t &vdof_data) const
284  {
285  const int NC = vdof_layout_t::dim_2;
286  const int NE = vdof_layout_t::dim_3;
287  MFEM_STATIC_ASSERT(FE::dofs == vdof_layout_t::dim_1,
288  "invalid number of dofs");
289  MFEM_ASSERT(first_comp + NC <= vl.NumComponents(),
290  "invalid number of components");
291  for (int k = 0; k < NC; k++)
292  {
293  for (int j = 0; j < NE; j++)
294  {
295  for (int i = 0; i < FE::dofs; i++)
296  {
297  Assign<AssignOp::Set>(
298  vdof_data[vdof_layout.ind(i,k,j)],
299  glob_vdof_data[vl.ind(ind.map(i,j), first_comp+k)]);
300  }
301  }
302  }
303  }
304 
305  // Assemble a static number of consecutive components; vdof_layout is
306  // (dofs x NC x NE), where NC is the number of components to add. It is
307  // assumed that: first_comp + NC <= vl.NumComponents().
308  template <typename vdof_layout_t, typename vdof_data_t,
309  typename vec_layout_t, typename glob_vdof_data_t>
310  inline MFEM_ALWAYS_INLINE
311  void AssembleComponents(int first_comp,
312  const vdof_layout_t &vdof_layout,
313  const vdof_data_t &vdof_data,
314  const vec_layout_t &vl,
315  glob_vdof_data_t &glob_vdof_data) const
316  {
317  const int NC = vdof_layout_t::dim_2;
318  const int NE = vdof_layout_t::dim_3;
319  MFEM_STATIC_ASSERT(FE::dofs == vdof_layout_t::dim_1,
320  "invalid number of dofs");
321  MFEM_ASSERT(first_comp + NC <= vl.NumComponents(),
322  "invalid number of components");
323  for (int k = 0; k < NC; k++)
324  {
325  for (int j = 0; j < NE; j++)
326  {
327  for (int i = 0; i < FE::dofs; i++)
328  {
329  Assign<AssignOp::Add>(
330  glob_vdof_data[vl.ind(ind.map(i,j), first_comp+k)],
331  vdof_data[vdof_layout.ind(i,k,j)]);
332  }
333  }
334  }
335  }
336 
338  SparseMatrix &M) const
339  {
340  MFEM_FLOPS_ADD(FE::dofs*FE::dofs);
341  for (int i = 0; i < FE::dofs; i++)
342  {
343  M.SetColPtr(ind.map(i,0));
344  for (int j = 0; j < FE::dofs; j++)
345  {
346  M._Add_(ind.map(j,0), m(i,j));
347  }
348  M.ClearColPtr();
349  }
350  }
351 
352  template <typename vec_layout_t>
353  void AssembleBlock(int block_i, int block_j, const vec_layout_t &vl,
355  SparseMatrix &M) const
356  {
357  MFEM_FLOPS_ADD(FE::dofs*FE::dofs);
358  for (int i = 0; i < FE::dofs; i++)
359  {
360  M.SetColPtr(vl.ind(ind.map(i,0), block_i));
361  for (int j = 0; j < FE::dofs; j++)
362  {
363  M._Add_(vl.ind(ind.map(j,0), block_j), m(i,j));
364  }
365  M.ClearColPtr();
366  }
367  }
368 };
369 
370 // H1 Finite Element Space
371 
372 template <typename FE>
374  : public TFiniteElementSpace_simple<FE,ElementDofIndexer<FE> >
375 {
376 public:
377  typedef FE FE_type;
379 
380  H1_FiniteElementSpace(const FE &fe, const FiniteElementSpace &fes)
381  : base_class(fe, fes)
382  { }
383 
384  // default copy constructor
385 
386  static bool Matches(const FiniteElementSpace &fes)
387  {
388  const FiniteElementCollection *fec = fes.FEColl();
389  const H1_FECollection *h1_fec =
390  dynamic_cast<const H1_FECollection *>(fec);
391  if (!h1_fec) { return false; }
393  if (fe->GetOrder() != FE_type::degree) { return false; }
394  return true;
395  }
396 
397  template <typename vec_layout_t>
398  static bool VectorMatches(const FiniteElementSpace &fes)
399  {
400  return Matches(fes) && vec_layout_t::Matches(fes);
401  }
402 };
403 
404 
405 // Simple index type for DG spaces, where the map method is given by:
406 // glob_dof_idx = loc_dof_idx + elem_idx * num_dofs.
407 template <typename FE>
409 {
410 protected:
411  int offset;
412 
413 public:
414  typedef FE FE_type;
415 
416  DGIndexer(const FE &fe, const FiniteElementSpace &fes)
417  {
418  MFEM_ASSERT(fes.GetNDofs() == fes.GetNE() * FE::dofs,
419  "the FE space is not compatible with this FE!");
420  offset = 0;
421  }
422 
423  // default copy constructor
424 
425  inline MFEM_ALWAYS_INLINE
426  void SetElement(int elem_idx)
427  {
428  offset = FE::dofs * elem_idx;
429  }
430 
431  inline MFEM_ALWAYS_INLINE
432  int map(int loc_dof_idx, int elem_offset) const
433  {
434  return offset + loc_dof_idx + elem_offset * FE::dofs;
435  }
436 };
437 
438 
439 // L2 Finite Element Space
440 
441 template <typename FE>
443  : public TFiniteElementSpace_simple<FE,DGIndexer<FE> >
444 {
445 public:
446  typedef FE FE_type;
448 
449  L2_FiniteElementSpace(const FE &fe, const FiniteElementSpace &fes)
450  : base_class(fe, fes)
451  { }
452 
453  // default copy constructor
454 
455  static bool Matches(const FiniteElementSpace &fes)
456  {
457  const FiniteElementCollection *fec = fes.FEColl();
458  const L2_FECollection *l2_fec =
459  dynamic_cast<const L2_FECollection *>(fec);
460  if (!l2_fec) { return false; }
462  if (fe->GetOrder() != FE_type::degree) { return false; }
463  return true;
464  }
465 
466  template <typename vec_layout_t>
467  static bool VectorMatches(const FiniteElementSpace &fes)
468  {
469  return Matches(fes) && vec_layout_t::Matches(fes);
470  }
471 };
472 
473 } // namespace mfem
474 
475 #endif // MFEM_TEMPLATE_FESPACE
Abstract class for Finite Elements.
Definition: fe.hpp:46
H1_FiniteElementSpace(const FE &fe, const FiniteElementSpace &fes)
Definition: tfespace.hpp:380
void _Add_(const int col, const double a)
Definition: sparsemat.hpp:297
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:161
int * GetJ()
Definition: table.hpp:108
MFEM_ALWAYS_INLINE void VectorExtract(const vec_layout_t &vl, const glob_vdof_data_t &glob_vdof_data, const vdof_layout_t &vdof_layout, vdof_data_t &vdof_data) const
Definition: tfespace.hpp:226
const Geometry::Type geom
MFEM_ALWAYS_INLINE int map(int loc_dof_idx, int elem_offset) const
Definition: tfespace.hpp:432
T * GetData()
Returns the data.
Definition: array.hpp:91
int GetOrder() const
Returns the order of the finite element.
Definition: fe.hpp:124
static bool VectorMatches(const FiniteElementSpace &fes)
Definition: tfespace.hpp:467
MFEM_ALWAYS_INLINE void ExtractComponents(int first_comp, const vec_layout_t &vl, const glob_vdof_data_t &glob_vdof_data, const vdof_layout_t &vdof_layout, vdof_data_t &vdof_data) const
Definition: tfespace.hpp:279
int Size_of_connections() const
Definition: table.hpp:92
const int * loc_dof_list
Definition: tfespace.hpp:43
MFEM_ALWAYS_INLINE void SetElement(int elem_idx)
Definition: tfespace.hpp:96
void Assemble(const TMatrix< FE::dofs, FE::dofs, double > &m, SparseMatrix &M) const
Definition: tfespace.hpp:337
MFEM_ALWAYS_INLINE void Extract(const glob_dof_data_t &glob_dof_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Definition: tfespace.hpp:156
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:184
void ClearColPtr() const
Definition: sparsemat.hpp:500
MFEM_ALWAYS_INLINE void VectorExtract(const vec_layout_t &vl, const glob_vdof_data_t &glob_vdof_data, const vdof_layout_t &vdof_layout, vdof_data_t &vdof_data) const
Definition: tfespace.hpp:200
L2_FiniteElementSpace(const FE &fe, const FiniteElementSpace &fes)
Definition: tfespace.hpp:449
Data type sparse matrix.
Definition: sparsemat.hpp:38
ElementDofIndexer(const FE &fe, const FiniteElementSpace &fes)
Definition: tfespace.hpp:49
MFEM_ALWAYS_INLINE void Assemble(const dof_layout_t &dof_layout, const dof_data_t &dof_data, glob_dof_data_t &glob_dof_data) const
Definition: tfespace.hpp:168
MFEM_ALWAYS_INLINE void Extract(const glob_dof_data_t &glob_dof_data, const dof_layout_t &dof_layout, dof_data_t &dof_data) const
Definition: tfespace.hpp:136
static bool VectorMatches(const FiniteElementSpace &fes)
Definition: tfespace.hpp:398
MFEM_ALWAYS_INLINE void SetElement(int elem_idx)
Definition: tfespace.hpp:426
MFEM_ALWAYS_INLINE void VectorAssemble(const vdof_layout_t &vdof_layout, const vdof_data_t &vdof_data, const vec_layout_t &vl, glob_vdof_data_t &glob_vdof_data) const
Definition: tfespace.hpp:265
const int * el_dof_list
Definition: tfespace.hpp:43
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:86
static bool Matches(const FiniteElementSpace &fes)
Definition: tfespace.hpp:455
static bool Matches(const FiniteElementSpace &fes)
Definition: tfespace.hpp:386
MFEM_ALWAYS_INLINE void AssembleComponents(int first_comp, const vdof_layout_t &vdof_layout, const vdof_data_t &vdof_data, const vec_layout_t &vl, glob_vdof_data_t &glob_vdof_data) const
Definition: tfespace.hpp:311
void SetColPtr(const int row) const
Definition: sparsemat.hpp:465
TFiniteElementSpace_simple< FE, ElementDofIndexer< FE > > base_class
Definition: tfespace.hpp:378
DGIndexer(const FE &fe, const FiniteElementSpace &fes)
Definition: tfespace.hpp:416
TFiniteElementSpace_simple(const FE &fe, const FiniteElementSpace &fes)
Definition: tfespace.hpp:122
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.hpp:160
MFEM_ALWAYS_INLINE ~ElementDofIndexer()
Definition: tfespace.hpp:93
virtual const FiniteElement * FiniteElementForGeometry(int GeomType) const
Definition: fe_coll.hpp:211
TFiniteElementSpace_simple< FE, DGIndexer< FE > > base_class
Definition: tfespace.hpp:447
void AssembleBlock(int block_i, int block_j, const vec_layout_t &vl, const TMatrix< FE::dofs, FE::dofs, double > &m, SparseMatrix &M) const
Definition: tfespace.hpp:353
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:177
MFEM_ALWAYS_INLINE ElementDofIndexer(const ElementDofIndexer &orig)
Definition: tfespace.hpp:86
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:146
const Table & GetElementToDofTable() const
Definition: fespace.hpp:288
int RowSize(int i) const
Definition: table.hpp:102
MFEM_ALWAYS_INLINE int map(int loc_dof_idx, int elem_offset) const
Definition: tfespace.hpp:102
MFEM_ALWAYS_INLINE void Assemble(const dof_layout_t &dof_layout, const dof_data_t &dof_data, glob_dof_data_t &glob_dof_data) const
Definition: tfespace.hpp:188
MFEM_ALWAYS_INLINE void VectorAssemble(const vdof_layout_t &vdof_layout, const vdof_data_t &vdof_data, const vec_layout_t &vl, glob_vdof_data_t &glob_vdof_data) const
Definition: tfespace.hpp:239
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:195