MFEM  v4.1.0 Finite element discretization library
restriction.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
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 #include "restriction.hpp"
13 #include "gridfunc.hpp"
14 #include "fespace.hpp"
15 #include "../general/forall.hpp"
16
17 namespace mfem
18 {
19
21  : ne(fes.GetNE()),
22  vdim(fes.GetVDim()),
23  byvdim(fes.GetOrdering() == Ordering::byVDIM),
24  ndof(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
25  ndofs(fes.GetNDofs())
26 {
27  height = vdim*ne*ndof;
28  width = vdim*ne*ndof;
29 }
30
31 void L2ElementRestriction::Mult(const Vector &x, Vector &y) const
32 {
33  const int nd = ndof;
34  const int vd = vdim;
35  const bool t = byvdim;
36  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
37  auto d_y = Reshape(y.Write(), nd, vd, ne);
38  MFEM_FORALL(i, ndofs,
39  {
40  const int idx = i;
41  const int dof = idx % nd;
42  const int e = idx / nd;
43  for (int c = 0; c < vd; ++c)
44  {
45  d_y(dof, c, e) = d_x(t?c:idx, t?idx:c);
46  }
47  });
48 }
49
51 {
52  const int nd = ndof;
53  const int vd = vdim;
54  const bool t = byvdim;
55  auto d_x = Reshape(x.Read(), nd, vd, ne);
56  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
57  MFEM_FORALL(i, ndofs,
58  {
59  const int idx = i;
60  const int dof = idx % nd;
61  const int e = idx / nd;
62  for (int c = 0; c < vd; ++c)
63  {
64  d_y(t?c:idx,t?idx:c) = d_x(dof, c, e);
65  }
66  });
67 }
68
70  ElementDofOrdering e_ordering)
71  : fes(f),
72  ne(fes.GetNE()),
73  vdim(fes.GetVDim()),
74  byvdim(fes.GetOrdering() == Ordering::byVDIM),
75  ndofs(fes.GetNDofs()),
76  dof(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
77  nedofs(ne*dof),
78  offsets(ndofs+1),
79  indices(ne*dof),
80  gatherMap(ne*dof)
81 {
82  // Assuming all finite elements are the same.
83  height = vdim*ne*dof;
84  width = fes.GetVSize();
85  const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
86  const int *dof_map = NULL;
87  if (dof_reorder && ne > 0)
88  {
89  for (int e = 0; e < ne; ++e)
90  {
91  const FiniteElement *fe = fes.GetFE(e);
92  const TensorBasisElement* el =
93  dynamic_cast<const TensorBasisElement*>(fe);
94  if (el) { continue; }
95  mfem_error("Finite element not suitable for lexicographic ordering");
96  }
97  const FiniteElement *fe = fes.GetFE(0);
98  const TensorBasisElement* el =
99  dynamic_cast<const TensorBasisElement*>(fe);
100  const Array<int> &fe_dof_map = el->GetDofMap();
101  MFEM_VERIFY(fe_dof_map.Size() > 0, "invalid dof map");
102  dof_map = fe_dof_map.GetData();
103  }
104  const Table& e2dTable = fes.GetElementToDofTable();
105  const int* elementMap = e2dTable.GetJ();
106  // We will be keeping a count of how many local nodes point to its global dof
107  for (int i = 0; i <= ndofs; ++i)
108  {
109  offsets[i] = 0;
110  }
111  for (int e = 0; e < ne; ++e)
112  {
113  for (int d = 0; d < dof; ++d)
114  {
115  const int sgid = elementMap[dof*e + d]; // signed
116  const int gid = (sgid >= 0) ? sgid : -1 - sgid;
117  ++offsets[gid + 1];
118  }
119  }
120  // Aggregate to find offsets for each global dof
121  for (int i = 1; i <= ndofs; ++i)
122  {
123  offsets[i] += offsets[i - 1];
124  }
125  // For each global dof, fill in all local nodes that point to it
126  for (int e = 0; e < ne; ++e)
127  {
128  for (int d = 0; d < dof; ++d)
129  {
130  const int sdid = dof_reorder ? dof_map[d] : 0; // signed
131  const int did = (!dof_reorder)?d:(sdid >= 0 ? sdid : -1-sdid);
132  const int sgid = elementMap[dof*e + did]; // signed
133  const int gid = (sgid >= 0) ? sgid : -1-sgid;
134  const int lid = dof*e + d;
135  const bool plus = (sgid >= 0 && sdid >= 0) || (sgid < 0 && sdid < 0);
136  gatherMap[lid] = plus ? gid : -1-gid;
137  indices[offsets[gid]++] = plus ? lid : -1-lid;
138  }
139  }
140  // We shifted the offsets vector by 1 by using it as a counter.
141  // Now we shift it back.
142  for (int i = ndofs; i > 0; --i)
143  {
144  offsets[i] = offsets[i - 1];
145  }
146  offsets[0] = 0;
147 }
148
149 void ElementRestriction::Mult(const Vector& x, Vector& y) const
150 {
151  // Assumes all elements have the same number of dofs
152  const int nd = dof;
153  const int vd = vdim;
154  const bool t = byvdim;
155  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
156  auto d_y = Reshape(y.Write(), nd, vd, ne);
158  MFEM_FORALL(i, dof*ne,
159  {
160  const int gid = d_gatherMap[i];
161  const bool plus = gid >= 0;
162  const int j = plus ? gid : -1-gid;
163  for (int c = 0; c < vd; ++c)
164  {
165  const double dofValue = d_x(t?c:j, t?j:c);
166  d_y(i % nd, c, i / nd) = plus ? dofValue : -dofValue;
167  }
168  });
169 }
170
172 {
173  // Assumes all elements have the same number of dofs
174  const int nd = dof;
175  const int vd = vdim;
176  const bool t = byvdim;
179  auto d_x = Reshape(x.Read(), nd, vd, ne);
180  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
181  MFEM_FORALL(i, ndofs,
182  {
183  const int offset = d_offsets[i];
184  const int nextOffset = d_offsets[i + 1];
185  for (int c = 0; c < vd; ++c)
186  {
187  double dofValue = 0;
188  for (int j = offset; j < nextOffset; ++j)
189  {
190  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
191  dofValue += (d_indices[j] >= 0) ? d_x(idx_j % nd, c,
192  idx_j / nd) : -d_x(idx_j % nd, c, idx_j / nd);
193  }
194  d_y(t?c:i,t?i:c) = dofValue;
195  }
196  });
197 }
198
200 {
201  // Assumes all elements have the same number of dofs
202  const int nd = dof;
203  const int vd = vdim;
204  const bool t = byvdim;
207  auto d_x = Reshape(x.Read(), nd, vd, ne);
208  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
209  MFEM_FORALL(i, ndofs,
210  {
211  const int offset = d_offsets[i];
212  const int nextOffset = d_offsets[i + 1];
213  for (int c = 0; c < vd; ++c)
214  {
215  double dofValue = 0;
216  for (int j = offset; j < nextOffset; ++j)
217  {
218  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
219  dofValue += d_x(idx_j % nd, c, idx_j / nd);
220  }
221  d_y(t?c:i,t?i:c) = dofValue;
222  }
223  });
224 }
225
226 /// Return the face degrees of freedom returned in Lexicographic order.
227 void GetFaceDofs(const int dim, const int face_id,
228  const int dof1d, Array<int> &faceMap)
229 {
230  switch (dim)
231  {
232  case 1:
233  switch (face_id)
234  {
235  case 0: // WEST
236  faceMap[0] = 0;
237  break;
238  case 1: // EAST
239  faceMap[0] = dof1d-1;
240  break;
241  }
242  break;
243  case 2:
244  switch (face_id)
245  {
246  case 0: // SOUTH
247  for (int i = 0; i < dof1d; ++i)
248  {
249  faceMap[i] = i;
250  }
251  break;
252  case 1: // EAST
253  for (int i = 0; i < dof1d; ++i)
254  {
255  faceMap[i] = dof1d-1 + i*dof1d;
256  }
257  break;
258  case 2: // NORTH
259  for (int i = 0; i < dof1d; ++i)
260  {
261  faceMap[i] = (dof1d-1)*dof1d + i;
262  }
263  break;
264  case 3: // WEST
265  for (int i = 0; i < dof1d; ++i)
266  {
267  faceMap[i] = i*dof1d;
268  }
269  break;
270  }
271  break;
272  case 3:
273  switch (face_id)
274  {
275  case 0: // BOTTOM
276  for (int i = 0; i < dof1d; ++i)
277  {
278  for (int j = 0; j < dof1d; ++j)
279  {
280  faceMap[i+j*dof1d] = i + j*dof1d;
281  }
282  }
283  break;
284  case 1: // SOUTH
285  for (int i = 0; i < dof1d; ++i)
286  {
287  for (int j = 0; j < dof1d; ++j)
288  {
289  faceMap[i+j*dof1d] = i + j*dof1d*dof1d;
290  }
291  }
292  break;
293  case 2: // EAST
294  for (int i = 0; i < dof1d; ++i)
295  {
296  for (int j = 0; j < dof1d; ++j)
297  {
298  faceMap[i+j*dof1d] = dof1d-1 + i*dof1d + j*dof1d*dof1d;
299  }
300  }
301  break;
302  case 3: // NORTH
303  for (int i = 0; i < dof1d; ++i)
304  {
305  for (int j = 0; j < dof1d; ++j)
306  {
307  faceMap[i+j*dof1d] = (dof1d-1)*dof1d + i + j*dof1d*dof1d;
308  }
309  }
310  break;
311  case 4: // WEST
312  for (int i = 0; i < dof1d; ++i)
313  {
314  for (int j = 0; j < dof1d; ++j)
315  {
316  faceMap[i+j*dof1d] = i*dof1d + j*dof1d*dof1d;
317  }
318  }
319  break;
320  case 5: // TOP
321  for (int i = 0; i < dof1d; ++i)
322  {
323  for (int j = 0; j < dof1d; ++j)
324  {
325  faceMap[i+j*dof1d] = (dof1d-1)*dof1d*dof1d + i + j*dof1d;
326  }
327  }
328  break;
329  }
330  break;
331  }
332 }
333
335  const ElementDofOrdering e_ordering,
336  const FaceType type)
337  : fes(fes),
338  nf(fes.GetNFbyType(type)),
339  vdim(fes.GetVDim()),
340  byvdim(fes.GetOrdering() == Ordering::byVDIM),
341  ndofs(fes.GetNDofs()),
342  dof(nf > 0 ? fes.GetFaceElement(0)->GetDof() : 0),
343  nfdofs(nf*dof),
344  scatter_indices(nf*dof),
345  offsets(ndofs+1),
346  gather_indices(nf*dof)
347 {
348  if (nf==0) { return; }
349  // If fespace == H1
350  const FiniteElement *fe = fes.GetFE(0);
351  const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
352  MFEM_VERIFY(tfe != NULL &&
355  "Only Gauss-Lobatto and Bernstein basis are supported in "
356  "H1FaceRestriction.");
357  MFEM_VERIFY(fes.GetMesh()->Conforming(),
358  "Non-conforming meshes not yet supported with partial assembly.");
359  // Assuming all finite elements are using Gauss-Lobatto.
360  height = vdim*nf*dof;
361  width = fes.GetVSize();
362  const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
363  if (dof_reorder && nf > 0)
364  {
365  for (int f = 0; f < fes.GetNF(); ++f)
366  {
367  const FiniteElement *fe = fes.GetFaceElement(f);
368  const TensorBasisElement* el =
369  dynamic_cast<const TensorBasisElement*>(fe);
370  if (el) { continue; }
371  mfem_error("Finite element not suitable for lexicographic ordering");
372  }
373  const FiniteElement *fe = fes.GetFaceElement(0);
374  const TensorBasisElement* el =
375  dynamic_cast<const TensorBasisElement*>(fe);
376  const Array<int> &fe_dof_map = el->GetDofMap();
377  MFEM_VERIFY(fe_dof_map.Size() > 0, "invalid dof map");
378  }
379  const TensorBasisElement* el =
380  dynamic_cast<const TensorBasisElement*>(fe);
381  const int *dof_map = el->GetDofMap().GetData();
382  const Table& e2dTable = fes.GetElementToDofTable();
383  const int* elementMap = e2dTable.GetJ();
384  Array<int> faceMap(dof);
385  int e1, e2;
386  int inf1, inf2;
387  int face_id;
388  int orientation;
389  const int dof1d = fes.GetFE(0)->GetOrder()+1;
390  const int elem_dofs = fes.GetFE(0)->GetDof();
391  const int dim = fes.GetMesh()->SpaceDimension();
392  // Computation of scatter_indices
393  int f_ind = 0;
394  for (int f = 0; f < fes.GetNF(); ++f)
395  {
396  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
397  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
398  orientation = inf1 % 64;
399  face_id = inf1 / 64;
400  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
401  (type==FaceType::Boundary && e2<0 && inf2<0) )
402  {
403  // Assumes Gauss-Lobatto basis
404  if (dof_reorder)
405  {
406  if (orientation != 0)
407  {
408  mfem_error("FaceRestriction used on degenerated mesh.");
409  }
410  GetFaceDofs(dim, face_id, dof1d, faceMap); // Only for hex
411  }
412  else
413  {
414  mfem_error("FaceRestriction not yet implemented for this type of "
415  "element.");
416  // TODO Something with GetFaceDofs?
417  }
418  for (int d = 0; d < dof; ++d)
419  {
420  const int face_dof = faceMap[d];
421  const int did = (!dof_reorder)?face_dof:dof_map[face_dof];
422  const int gid = elementMap[e1*elem_dofs + did];
423  const int lid = dof*f_ind + d;
424  scatter_indices[lid] = gid;
425  }
426  f_ind++;
427  }
428  }
429  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
430  // Computation of gather_indices
431  for (int i = 0; i <= ndofs; ++i)
432  {
433  offsets[i] = 0;
434  }
435  f_ind = 0;
436  for (int f = 0; f < fes.GetNF(); ++f)
437  {
438  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
439  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
440  orientation = inf1 % 64;
441  face_id = inf1 / 64;
442  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
443  (type==FaceType::Boundary && e2<0 && inf2<0) )
444  {
445  GetFaceDofs(dim, face_id, dof1d, faceMap);
446  for (int d = 0; d < dof; ++d)
447  {
448  const int face_dof = faceMap[d];
449  const int did = (!dof_reorder)?face_dof:dof_map[face_dof];
450  const int gid = elementMap[e1*elem_dofs + did];
451  ++offsets[gid + 1];
452  }
453  f_ind++;
454  }
455  }
456  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
457  for (int i = 1; i <= ndofs; ++i)
458  {
459  offsets[i] += offsets[i - 1];
460  }
461  f_ind = 0;
462  for (int f = 0; f < fes.GetNF(); ++f)
463  {
464  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
465  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
466  orientation = inf1 % 64;
467  face_id = inf1 / 64;
468  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
469  (type==FaceType::Boundary && e2<0 && inf2<0) )
470  {
471  GetFaceDofs(dim, face_id, dof1d, faceMap);
472  for (int d = 0; d < dof; ++d)
473  {
474  const int face_dof = faceMap[d];
475  const int did = (!dof_reorder)?face_dof:dof_map[face_dof];
476  const int gid = elementMap[e1*elem_dofs + did];
477  const int lid = dof*f_ind + d;
478  gather_indices[offsets[gid]++] = lid;
479  }
480  f_ind++;
481  }
482  }
483  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
484  for (int i = ndofs; i > 0; --i)
485  {
486  offsets[i] = offsets[i - 1];
487  }
488  offsets[0] = 0;
489 }
490
491 void H1FaceRestriction::Mult(const Vector& x, Vector& y) const
492 {
493  // Assumes all elements have the same number of dofs
494  const int nd = dof;
495  const int vd = vdim;
496  const bool t = byvdim;
498  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
499  auto d_y = Reshape(y.Write(), nd, vd, nf);
500  MFEM_FORALL(i, nfdofs,
501  {
502  const int idx = d_indices[i];
503  const int dof = i % nd;
504  const int face = i / nd;
505  for (int c = 0; c < vd; ++c)
506  {
507  d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
508  }
509  });
510 }
511
513 {
514  // Assumes all elements have the same number of dofs
515  const int nd = dof;
516  const int vd = vdim;
517  const bool t = byvdim;
520  auto d_x = Reshape(x.Read(), nd, vd, nf);
521  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
522  MFEM_FORALL(i, ndofs,
523  {
524  const int offset = d_offsets[i];
525  const int nextOffset = d_offsets[i + 1];
526  for (int c = 0; c < vd; ++c)
527  {
528  double dofValue = 0;
529  for (int j = offset; j < nextOffset; ++j)
530  {
531  const int idx_j = d_indices[j];
532  dofValue += d_x(idx_j % nd, c, idx_j / nd);
533  }
534  d_y(t?c:i,t?i:c) += dofValue;
535  }
536  });
537 }
538
539 static int ToLexOrdering2D(const int face_id, const int size1d, const int i)
540 {
541  if (face_id==2 || face_id==3)
542  {
543  return size1d-1-i;
544  }
545  else
546  {
547  return i;
548  }
549 }
550
551 static int PermuteFace2D(const int face_id1, const int face_id2,
552  const int orientation,
553  const int size1d, const int index)
554 {
555  int new_index;
556  // Convert from lex ordering
557  if (face_id1==2 || face_id1==3)
558  {
559  new_index = size1d-1-index;
560  }
561  else
562  {
563  new_index = index;
564  }
565  // Permute based on face orientations
566  if (orientation==1)
567  {
568  new_index = size1d-1-new_index;
569  }
571 }
572
573 static int ToLexOrdering3D(const int face_id, const int size1d, const int i,
574  const int j)
575 {
576  if (face_id==2 || face_id==1 || face_id==5)
577  {
578  return i + j*size1d;
579  }
580  else if (face_id==3 || face_id==4)
581  {
582  return (size1d-1-i) + j*size1d;
583  }
584  else // face_id==0
585  {
586  return i + (size1d-1-j)*size1d;
587  }
588 }
589
590 static int PermuteFace3D(const int face_id1, const int face_id2,
591  const int orientation,
592  const int size1d, const int index)
593 {
594  int i=0, j=0, new_i=0, new_j=0;
595  i = index%size1d;
596  j = index/size1d;
597  // Convert from lex ordering
598  if (face_id1==3 || face_id1==4)
599  {
600  i = size1d-1-i;
601  }
602  else if (face_id1==0)
603  {
604  j = size1d-1-j;
605  }
606  // Permute based on face orientations
607  switch (orientation)
608  {
609  case 0:
610  new_i = i;
611  new_j = j;
612  break;
613  case 1:
614  new_i = j;
615  new_j = i;
616  break;
617  case 2:
618  new_i = j;
619  new_j = (size1d-1-i);
620  break;
621  case 3:
622  new_i = (size1d-1-i);
623  new_j = j;
624  break;
625  case 4:
626  new_i = (size1d-1-i);
627  new_j = (size1d-1-j);
628  break;
629  case 5:
630  new_i = (size1d-1-j);
631  new_j = (size1d-1-i);
632  break;
633  case 6:
634  new_i = (size1d-1-j);
635  new_j = i;
636  break;
637  case 7:
638  new_i = i;
639  new_j = (size1d-1-j);
640  break;
641  }
643 }
644
645 /// Permute dofs or quads on a face for e2 to match with the ordering of e1
646 int PermuteFaceL2(const int dim, const int face_id1,
647  const int face_id2, const int orientation,
648  const int size1d, const int index)
649 {
650  switch (dim)
651  {
652  case 1:
653  return 0;
654  case 2:
655  return PermuteFace2D(face_id1, face_id2, orientation, size1d, index);
656  case 3:
657  return PermuteFace3D(face_id1, face_id2, orientation, size1d, index);
658  default:
659  mfem_error("Unsupported dimension.");
660  return 0;
661  }
662 }
663
665  const ElementDofOrdering e_ordering,
666  const FaceType type,
667  const L2FaceValues m)
668  : fes(fes),
669  nf(fes.GetNFbyType(type)),
670  vdim(fes.GetVDim()),
671  byvdim(fes.GetOrdering() == Ordering::byVDIM),
672  ndofs(fes.GetNDofs()),
673  dof(nf > 0 ?
674  fes.GetTraceElement(0, fes.GetMesh()->GetFaceBaseGeometry(0))->GetDof()
675  : 0),
676  m(m),
677  nfdofs(nf*dof),
678  scatter_indices1(nf*dof),
679  scatter_indices2(m==L2FaceValues::DoubleValued?nf*dof:0),
680  offsets(ndofs+1),
681  gather_indices((m==L2FaceValues::DoubleValued? 2 : 1)*nf*dof)
682 {
683  // If fespace == L2
684  const FiniteElement *fe = fes.GetFE(0);
685  const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
686  MFEM_VERIFY(tfe != NULL &&
689  "Only Gauss-Lobatto and Bernstein basis are supported in "
690  "L2FaceRestriction.");
691  MFEM_VERIFY(fes.GetMesh()->Conforming(),
692  "Non-conforming meshes not yet supported with partial assembly.");
693  if (nf==0) { return; }
695  width = fes.GetVSize();
696  const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
697  if (!dof_reorder)
698  {
699  mfem_error("Non-Tensor L2FaceRestriction not yet implemented.");
700  }
701  if (dof_reorder && nf > 0)
702  {
703  for (int f = 0; f < fes.GetNF(); ++f)
704  {
705  const FiniteElement *fe =
706  fes.GetTraceElement(f, fes.GetMesh()->GetFaceBaseGeometry(f));
707  const TensorBasisElement* el =
708  dynamic_cast<const TensorBasisElement*>(fe);
709  if (el) { continue; }
710  mfem_error("Finite element not suitable for lexicographic ordering");
711  }
712  }
713  const Table& e2dTable = fes.GetElementToDofTable();
714  const int* elementMap = e2dTable.GetJ();
715  Array<int> faceMap1(dof), faceMap2(dof);
716  int e1, e2;
717  int inf1, inf2;
718  int face_id1, face_id2;
719  int orientation;
720  const int dof1d = fes.GetFE(0)->GetOrder()+1;
721  const int elem_dofs = fes.GetFE(0)->GetDof();
722  const int dim = fes.GetMesh()->SpaceDimension();
723  // Computation of scatter indices
724  int f_ind=0;
725  for (int f = 0; f < fes.GetNF(); ++f)
726  {
727  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
728  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
729  if (dof_reorder)
730  {
731  orientation = inf1 % 64;
732  face_id1 = inf1 / 64;
733  GetFaceDofs(dim, face_id1, dof1d, faceMap1); // Only for hex
734  orientation = inf2 % 64;
735  face_id2 = inf2 / 64;
736  GetFaceDofs(dim, face_id2, dof1d, faceMap2); // Only for hex
737  }
738  else
739  {
740  mfem_error("FaceRestriction not yet implemented for this type of "
741  "element.");
742  // TODO Something with GetFaceDofs?
743  }
744  if ((type==FaceType::Interior && e2>=0) ||
745  (type==FaceType::Boundary && e2<0))
746  {
747  for (int d = 0; d < dof; ++d)
748  {
749  const int face_dof = faceMap1[d];
750  const int did = face_dof;
751  const int gid = elementMap[e1*elem_dofs + did];
752  const int lid = dof*f_ind + d;
753  scatter_indices1[lid] = gid;
754  }
756  {
757  for (int d = 0; d < dof; ++d)
758  {
759  if (type==FaceType::Interior && e2>=0) // interior face
760  {
761  const int pd = PermuteFaceL2(dim, face_id1, face_id2,
762  orientation, dof1d, d);
763  const int face_dof = faceMap2[pd];
764  const int did = face_dof;
765  const int gid = elementMap[e2*elem_dofs + did];
766  const int lid = dof*f_ind + d;
767  scatter_indices2[lid] = gid;
768  }
769  else if (type==FaceType::Boundary && e2<0) // true boundary face
770  {
771  const int lid = dof*f_ind + d;
772  scatter_indices2[lid] = -1;
773  }
774  }
775  }
776  f_ind++;
777  }
778  }
779  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
780  // Computation of gather_indices
781  for (int i = 0; i <= ndofs; ++i)
782  {
783  offsets[i] = 0;
784  }
785  f_ind = 0;
786  for (int f = 0; f < fes.GetNF(); ++f)
787  {
788  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
789  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
790  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
791  (type==FaceType::Boundary && e2<0 && inf2<0) )
792  {
793  orientation = inf1 % 64;
794  face_id1 = inf1 / 64;
795  GetFaceDofs(dim, face_id1, dof1d, faceMap1);
796  orientation = inf2 % 64;
797  face_id2 = inf2 / 64;
798  GetFaceDofs(dim, face_id2, dof1d, faceMap2);
799  for (int d = 0; d < dof; ++d)
800  {
801  const int did = faceMap1[d];
802  const int gid = elementMap[e1*elem_dofs + did];
803  ++offsets[gid + 1];
804  }
806  {
807  for (int d = 0; d < dof; ++d)
808  {
809  if (type==FaceType::Interior && e2>=0) // interior face
810  {
811  const int pd = PermuteFaceL2(dim, face_id1, face_id2,
812  orientation, dof1d, d);
813  const int did = faceMap2[pd];
814  const int gid = elementMap[e2*elem_dofs + did];
815  ++offsets[gid + 1];
816  }
817  }
818  }
819  f_ind++;
820  }
821  }
822  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
823  for (int i = 1; i <= ndofs; ++i)
824  {
825  offsets[i] += offsets[i - 1];
826  }
827  f_ind = 0;
828  for (int f = 0; f < fes.GetNF(); ++f)
829  {
830  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
831  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
832  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
833  (type==FaceType::Boundary && e2<0 && inf2<0) )
834  {
835  orientation = inf1 % 64;
836  face_id1 = inf1 / 64;
837  GetFaceDofs(dim, face_id1, dof1d, faceMap1);
838  orientation = inf2 % 64;
839  face_id2 = inf2 / 64;
840  GetFaceDofs(dim, face_id2, dof1d, faceMap2);
841  for (int d = 0; d < dof; ++d)
842  {
843  const int did = faceMap1[d];
844  const int gid = elementMap[e1*elem_dofs + did];
845  const int lid = dof*f_ind + d;
846  // We don't shift lid to express that it's e1 of f
847  gather_indices[offsets[gid]++] = lid;
848  }
850  {
851  for (int d = 0; d < dof; ++d)
852  {
853  if (type==FaceType::Interior && e2>=0) // interior face
854  {
855  const int pd = PermuteFaceL2(dim, face_id1, face_id2,
856  orientation, dof1d, d);
857  const int did = faceMap2[pd];
858  const int gid = elementMap[e2*elem_dofs + did];
859  const int lid = dof*f_ind + d;
860  // We shift lid to express that it's e2 of f
861  gather_indices[offsets[gid]++] = nfdofs + lid;
862  }
863  }
864  }
865  f_ind++;
866  }
867  }
868  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
869  for (int i = ndofs; i > 0; --i)
870  {
871  offsets[i] = offsets[i - 1];
872  }
873  offsets[0] = 0;
874 }
875
876 void L2FaceRestriction::Mult(const Vector& x, Vector& y) const
877 {
878  // Assumes all elements have the same number of dofs
879  const int nd = dof;
880  const int vd = vdim;
881  const bool t = byvdim;
882
884  {
887  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
888  auto d_y = Reshape(y.Write(), nd, vd, 2, nf);
889  MFEM_FORALL(i, nfdofs,
890  {
891  const int dof = i % nd;
892  const int face = i / nd;
893  const int idx1 = d_indices1[i];
894  for (int c = 0; c < vd; ++c)
895  {
896  d_y(dof, c, 0, face) = d_x(t?c:idx1, t?idx1:c);
897  }
898  const int idx2 = d_indices2[i];
899  for (int c = 0; c < vd; ++c)
900  {
901  d_y(dof, c, 1, face) = idx2==-1 ? 0.0 : d_x(t?c:idx2, t?idx2:c);
902  }
903  });
904  }
905  else
906  {
908  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
909  auto d_y = Reshape(y.Write(), nd, vd, nf);
910  MFEM_FORALL(i, nfdofs,
911  {
912  const int dof = i % nd;
913  const int face = i / nd;
914  const int idx1 = d_indices1[i];
915  for (int c = 0; c < vd; ++c)
916  {
917  d_y(dof, c, face) = d_x(t?c:idx1, t?idx1:c);
918  }
919  });
920  }
921 }
922
924 {
925  // Assumes all elements have the same number of dofs
926  const int nd = dof;
927  const int vd = vdim;
928  const bool t = byvdim;
929  const int dofs = nfdofs;
932  auto d_x = Reshape(x.Read(), nd, vd, 2, nf);
933  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
934  MFEM_FORALL(i, ndofs,
935  {
936  const int offset = d_offsets[i];
937  const int nextOffset = d_offsets[i + 1];
938  for (int c = 0; c < vd; ++c)
939  {
940  double dofValue = 0;
941  for (int j = offset; j < nextOffset; ++j)
942  {
943  int idx_j = d_indices[j];
944  bool isE1 = idx_j < dofs;
945  idx_j = isE1 ? idx_j : idx_j - dofs;
946  dofValue += isE1 ?
947  d_x(idx_j % nd, c, 0, idx_j / nd)
948  :d_x(idx_j % nd, c, 1, idx_j / nd);
949  }
950  d_y(t?c:i,t?i:c) += dofValue;
951  }
952  });
953 }
954
955 int ToLexOrdering(const int dim, const int face_id, const int size1d,
956  const int index)
957 {
958  switch (dim)
959  {
960  case 1:
961  return 0;
962  case 2:
964  case 3:
966  default:
967  mfem_error("Unsupported dimension.");
968  return 0;
969  }
970 }
971
972 } // namespace mfem
Abstract class for Finite Elements.
Definition: fe.hpp:232
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
int Size() const
Logical size of the array.
Definition: array.hpp:124
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:381
int * GetJ()
Definition: table.hpp:114
L2FaceValues
Definition: restriction.hpp:26
void MultTransposeUnsigned(const Vector &x, Vector &y) const
Compute MultTranspose without applying signs based on DOF orientations.
L2FaceRestriction(const FiniteElementSpace &, const ElementDofOrdering, const FaceType, const L2FaceValues m=L2FaceValues::DoubleValued)
const FiniteElement * GetTraceElement(int i, Geometry::Type geom_type) const
Return the trace element from element &#39;i&#39; to the given &#39;geom_type&#39;.
Definition: fespace.cpp:1936
bool Conforming() const
Definition: mesh.hpp:1167
Array< int > gather_indices
Definition: restriction.hpp:87
T * GetData()
Returns the data.
Definition: array.hpp:98
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:318
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
H1FaceRestriction(const FiniteElementSpace &, const ElementDofOrdering, const FaceType)
Geometry::Type GetFaceBaseGeometry(int i) const
Definition: mesh.hpp:785
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
const FiniteElement * GetFaceElement(int i) const
Definition: fespace.cpp:1908
double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:341
const L2FaceValues m
Array< int > gather_indices
L2ElementRestriction(const FiniteElementSpace &)
Definition: restriction.cpp:20
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:134
int GetBasisType() const
Definition: fe.hpp:1842
FaceType
Definition: mesh.hpp:42
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: restriction.cpp:50
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:295
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
const T * Read(bool on_dev=true) const
Definition: array.hpp:273
ElementRestriction(const FiniteElementSpace &, ElementDofOrdering)
Definition: restriction.cpp:69
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
Definition: fespace.hpp:414
int SpaceDimension() const
Definition: mesh.hpp:745
int ToLexOrdering(const int dim, const int face_id, const int size1d, const int index)
const FiniteElementSpace & fes
Definition: restriction.hpp:34
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:974
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: restriction.cpp:31
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:980
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:314
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void GetFaceDofs(const int dim, const int face_id, const int dof1d, Array< int > &faceMap)
Return the face degrees of freedom returned in Lexicographic order.
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition: fespace.hpp:28
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
const Array< int > & GetDofMap() const
Get an Array&lt;int&gt; that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe.hpp:1850
Array< int > scatter_indices1
const double * Read(bool on_dev=true) const
Definition: vector.hpp:333
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:65
int dim
Definition: ex24.cpp:43
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
int index(int i, int j, int nx, int ny)
Definition: life.cpp:241
Lexicographic ordering for tensor-product FiniteElements.
Vector data type.
Definition: vector.hpp:48
int PermuteFaceL2(const int dim, const int face_id1, const int face_id2, const int orientation, const int size1d, const int index)
Permute dofs or quads on a face for e2 to match with the ordering of e1.
const Table & GetElementToDofTable() const
Definition: fespace.hpp:524
Bernstein polynomials.
Definition: fe.hpp:35
Array< int > scatter_indices
Definition: restriction.hpp:85
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
Array< int > scatter_indices2
VDIM x NQPT x NE.
void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...