MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
restriction.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 #include "restriction.hpp"
13 #include "gridfunc.hpp"
14 #include "fespace.hpp"
15 #include "../general/forall.hpp"
16 #include <climits>
17 
18 namespace mfem
19 {
20 
22  ElementDofOrdering e_ordering)
23  : fes(f),
24  ne(fes.GetNE()),
25  vdim(fes.GetVDim()),
26  byvdim(fes.GetOrdering() == Ordering::byVDIM),
27  ndofs(fes.GetNDofs()),
28  dof(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
29  nedofs(ne*dof),
30  offsets(ndofs+1),
31  indices(ne*dof),
32  gatherMap(ne*dof)
33 {
34  // Assuming all finite elements are the same.
35  height = vdim*ne*dof;
36  width = fes.GetVSize();
37  const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
38  const int *dof_map = NULL;
39  if (dof_reorder && ne > 0)
40  {
41  for (int e = 0; e < ne; ++e)
42  {
43  const FiniteElement *fe = fes.GetFE(e);
44  const TensorBasisElement* el =
45  dynamic_cast<const TensorBasisElement*>(fe);
46  if (el) { continue; }
47  mfem_error("Finite element not suitable for lexicographic ordering");
48  }
49  const FiniteElement *fe = fes.GetFE(0);
50  const TensorBasisElement* el =
51  dynamic_cast<const TensorBasisElement*>(fe);
52  const Array<int> &fe_dof_map = el->GetDofMap();
53  MFEM_VERIFY(fe_dof_map.Size() > 0, "invalid dof map");
54  dof_map = fe_dof_map.GetData();
55  }
56  const Table& e2dTable = fes.GetElementToDofTable();
57  const int* elementMap = e2dTable.GetJ();
58  // We will be keeping a count of how many local nodes point to its global dof
59  for (int i = 0; i <= ndofs; ++i)
60  {
61  offsets[i] = 0;
62  }
63  for (int e = 0; e < ne; ++e)
64  {
65  for (int d = 0; d < dof; ++d)
66  {
67  const int sgid = elementMap[dof*e + d]; // signed
68  const int gid = (sgid >= 0) ? sgid : -1 - sgid;
69  ++offsets[gid + 1];
70  }
71  }
72  // Aggregate to find offsets for each global dof
73  for (int i = 1; i <= ndofs; ++i)
74  {
75  offsets[i] += offsets[i - 1];
76  }
77  // For each global dof, fill in all local nodes that point to it
78  for (int e = 0; e < ne; ++e)
79  {
80  for (int d = 0; d < dof; ++d)
81  {
82  const int sdid = dof_reorder ? dof_map[d] : 0; // signed
83  const int did = (!dof_reorder)?d:(sdid >= 0 ? sdid : -1-sdid);
84  const int sgid = elementMap[dof*e + did]; // signed
85  const int gid = (sgid >= 0) ? sgid : -1-sgid;
86  const int lid = dof*e + d;
87  const bool plus = (sgid >= 0 && sdid >= 0) || (sgid < 0 && sdid < 0);
88  gatherMap[lid] = plus ? gid : -1-gid;
89  indices[offsets[gid]++] = plus ? lid : -1-lid;
90  }
91  }
92  // We shifted the offsets vector by 1 by using it as a counter.
93  // Now we shift it back.
94  for (int i = ndofs; i > 0; --i)
95  {
96  offsets[i] = offsets[i - 1];
97  }
98  offsets[0] = 0;
99 }
100 
101 void ElementRestriction::Mult(const Vector& x, Vector& y) const
102 {
103  // Assumes all elements have the same number of dofs
104  const int nd = dof;
105  const int vd = vdim;
106  const bool t = byvdim;
107  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
108  auto d_y = Reshape(y.Write(), nd, vd, ne);
109  auto d_gatherMap = gatherMap.Read();
110  MFEM_FORALL(i, dof*ne,
111  {
112  const int gid = d_gatherMap[i];
113  const bool plus = gid >= 0;
114  const int j = plus ? gid : -1-gid;
115  for (int c = 0; c < vd; ++c)
116  {
117  const double dofValue = d_x(t?c:j, t?j:c);
118  d_y(i % nd, c, i / nd) = plus ? dofValue : -dofValue;
119  }
120  });
121 }
122 
124 {
125  // Assumes all elements have the same number of dofs
126  const int nd = dof;
127  const int vd = vdim;
128  const bool t = byvdim;
129  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
130  auto d_y = Reshape(y.Write(), nd, vd, ne);
131  auto d_gatherMap = gatherMap.Read();
132 
133  MFEM_FORALL(i, dof*ne,
134  {
135  const int gid = d_gatherMap[i];
136  const int j = gid >= 0 ? gid : -1-gid;
137  for (int c = 0; c < vd; ++c)
138  {
139  d_y(i % nd, c, i / nd) = d_x(t?c:j, t?j:c);
140  }
141  });
142 }
143 
145 {
146  // Assumes all elements have the same number of dofs
147  const int nd = dof;
148  const int vd = vdim;
149  const bool t = byvdim;
150  auto d_offsets = offsets.Read();
151  auto d_indices = indices.Read();
152  auto d_x = Reshape(x.Read(), nd, vd, ne);
153  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
154  MFEM_FORALL(i, ndofs,
155  {
156  const int offset = d_offsets[i];
157  const int nextOffset = d_offsets[i + 1];
158  for (int c = 0; c < vd; ++c)
159  {
160  double dofValue = 0;
161  for (int j = offset; j < nextOffset; ++j)
162  {
163  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
164  dofValue += (d_indices[j] >= 0) ? d_x(idx_j % nd, c,
165  idx_j / nd) : -d_x(idx_j % nd, c, idx_j / nd);
166  }
167  d_y(t?c:i,t?i:c) = dofValue;
168  }
169  });
170 }
171 
173 {
174  // Assumes all elements have the same number of dofs
175  const int nd = dof;
176  const int vd = vdim;
177  const bool t = byvdim;
178  auto d_offsets = offsets.Read();
179  auto d_indices = indices.Read();
180  auto d_x = Reshape(x.Read(), nd, vd, ne);
181  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
182  MFEM_FORALL(i, ndofs,
183  {
184  const int offset = d_offsets[i];
185  const int nextOffset = d_offsets[i + 1];
186  for (int c = 0; c < vd; ++c)
187  {
188  double dofValue = 0;
189  for (int j = offset; j < nextOffset; ++j)
190  {
191  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
192  dofValue += 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;
205  auto d_offsets = offsets.Read();
206  auto d_indices = indices.Read();
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 nextOffset = d_offsets[i + 1];
212  for (int c = 0; c < vd; ++c)
213  {
214  double dofValue = 0;
215  const int j = nextOffset - 1;
216  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
217  dofValue = (d_indices[j] >= 0) ? d_x(idx_j % nd, c, idx_j / nd) :
218  -d_x(idx_j % nd, c, idx_j / nd);
219  d_y(t?c:i,t?i:c) = dofValue;
220  }
221  });
222 }
223 
225 {
226  // Assumes all elements have the same number of dofs
227  const int nd = dof;
228  const int vd = vdim;
229  const bool t = byvdim;
230 
231  Array<char> processed(vd * ndofs);
232  processed = 0;
233 
234  auto d_offsets = offsets.HostRead();
235  auto d_indices = indices.HostRead();
236  auto d_x = Reshape(processed.HostReadWrite(), t?vd:ndofs, t?ndofs:vd);
237  auto d_y = Reshape(y.HostWrite(), nd, vd, ne);
238  for (int i = 0; i < ndofs; ++i)
239  {
240  const int offset = d_offsets[i];
241  const int nextOffset = d_offsets[i+1];
242  for (int c = 0; c < vd; ++c)
243  {
244  for (int j = offset; j < nextOffset; ++j)
245  {
246  const int idx_j = d_indices[j];
247  if (d_x(t?c:i,t?i:c))
248  {
249  d_y(idx_j % nd, c, idx_j / nd) = 0.0;
250  }
251  else
252  {
253  d_y(idx_j % nd, c, idx_j / nd) = 1.0;
254  d_x(t?c:i,t?i:c) = 1;
255  }
256  }
257  }
258  }
259 }
260 
262  SparseMatrix &mat) const
263 {
264  mat.GetMemoryI().New(mat.Height()+1, mat.GetMemoryI().GetMemoryType());
265  const int nnz = FillI(mat);
266  mat.GetMemoryJ().New(nnz, mat.GetMemoryJ().GetMemoryType());
267  mat.GetMemoryData().New(nnz, mat.GetMemoryData().GetMemoryType());
268  FillJAndData(mat_ea, mat);
269 }
270 
271 static MFEM_HOST_DEVICE int GetMinElt(const int *my_elts, const int nbElts,
272  const int *nbr_elts, const int nbrNbElts)
273 {
274  // Find the minimal element index found in both my_elts[] and nbr_elts[]
275  int min_el = INT_MAX;
276  for (int i = 0; i < nbElts; i++)
277  {
278  const int e_i = my_elts[i];
279  if (e_i >= min_el) { continue; }
280  for (int j = 0; j < nbrNbElts; j++)
281  {
282  if (e_i==nbr_elts[j])
283  {
284  min_el = e_i; // we already know e_i < min_el
285  break;
286  }
287  }
288  }
289  return min_el;
290 }
291 
292 /** Returns the index where a non-zero entry should be added and increment the
293  number of non-zeros for the row i_L. */
294 static MFEM_HOST_DEVICE int GetAndIncrementNnzIndex(const int i_L, int* I)
295 {
296  int ind = AtomicAdd(I[i_L],1);
297  return ind;
298 }
299 
301 {
302  static constexpr int Max = MaxNbNbr;
303  const int all_dofs = ndofs;
304  const int vd = vdim;
305  const int elt_dofs = dof;
306  auto I = mat.ReadWriteI();
307  auto d_offsets = offsets.Read();
308  auto d_indices = indices.Read();
309  auto d_gatherMap = gatherMap.Read();
310  MFEM_FORALL(i_L, vd*all_dofs+1,
311  {
312  I[i_L] = 0;
313  });
314  MFEM_FORALL(e, ne,
315  {
316  for (int i = 0; i < elt_dofs; i++)
317  {
318  int i_elts[Max];
319  const int i_E = e*elt_dofs + i;
320  const int i_L = d_gatherMap[i_E];
321  const int i_offset = d_offsets[i_L];
322  const int i_nextOffset = d_offsets[i_L+1];
323  const int i_nbElts = i_nextOffset - i_offset;
324  for (int e_i = 0; e_i < i_nbElts; ++e_i)
325  {
326  const int i_E = d_indices[i_offset+e_i];
327  i_elts[e_i] = i_E/elt_dofs;
328  }
329  for (int j = 0; j < elt_dofs; j++)
330  {
331  const int j_E = e*elt_dofs + j;
332  const int j_L = d_gatherMap[j_E];
333  const int j_offset = d_offsets[j_L];
334  const int j_nextOffset = d_offsets[j_L+1];
335  const int j_nbElts = j_nextOffset - j_offset;
336  if (i_nbElts == 1 || j_nbElts == 1) // no assembly required
337  {
338  GetAndIncrementNnzIndex(i_L, I);
339  }
340  else // assembly required
341  {
342  int j_elts[Max];
343  for (int e_j = 0; e_j < j_nbElts; ++e_j)
344  {
345  const int j_E = d_indices[j_offset+e_j];
346  const int elt = j_E/elt_dofs;
347  j_elts[e_j] = elt;
348  }
349  int min_e = GetMinElt(i_elts, i_nbElts, j_elts, j_nbElts);
350  if (e == min_e) // add the nnz only once
351  {
352  GetAndIncrementNnzIndex(i_L, I);
353  }
354  }
355  }
356  }
357  });
358  // We need to sum the entries of I, we do it on CPU as it is very sequential.
359  auto h_I = mat.HostReadWriteI();
360  const int nTdofs = vd*all_dofs;
361  int sum = 0;
362  for (int i = 0; i < nTdofs; i++)
363  {
364  const int nnz = h_I[i];
365  h_I[i] = sum;
366  sum+=nnz;
367  }
368  h_I[nTdofs] = sum;
369  // We return the number of nnz
370  return h_I[nTdofs];
371 }
372 
374  SparseMatrix &mat) const
375 {
376  static constexpr int Max = MaxNbNbr;
377  const int all_dofs = ndofs;
378  const int vd = vdim;
379  const int elt_dofs = dof;
380  auto I = mat.ReadWriteI();
381  auto J = mat.WriteJ();
382  auto Data = mat.WriteData();
383  auto d_offsets = offsets.Read();
384  auto d_indices = indices.Read();
385  auto d_gatherMap = gatherMap.Read();
386  auto mat_ea = Reshape(ea_data.Read(), elt_dofs, elt_dofs, ne);
387  MFEM_FORALL(e, ne,
388  {
389  for (int i = 0; i < elt_dofs; i++)
390  {
391  int i_elts[Max];
392  int i_B[Max];
393  const int i_E = e*elt_dofs + i;
394  const int i_L = d_gatherMap[i_E];
395  const int i_offset = d_offsets[i_L];
396  const int i_nextOffset = d_offsets[i_L+1];
397  const int i_nbElts = i_nextOffset - i_offset;
398  for (int e_i = 0; e_i < i_nbElts; ++e_i)
399  {
400  const int i_E = d_indices[i_offset+e_i];
401  i_elts[e_i] = i_E/elt_dofs;
402  i_B[e_i] = i_E%elt_dofs;
403  }
404  for (int j = 0; j < elt_dofs; j++)
405  {
406  const int j_E = e*elt_dofs + j;
407  const int j_L = d_gatherMap[j_E];
408  const int j_offset = d_offsets[j_L];
409  const int j_nextOffset = d_offsets[j_L+1];
410  const int j_nbElts = j_nextOffset - j_offset;
411  if (i_nbElts == 1 || j_nbElts == 1) // no assembly required
412  {
413  const int nnz = GetAndIncrementNnzIndex(i_L, I);
414  J[nnz] = j_L;
415  Data[nnz] = mat_ea(j,i,e);
416  }
417  else // assembly required
418  {
419  int j_elts[Max];
420  int j_B[Max];
421  for (int e_j = 0; e_j < j_nbElts; ++e_j)
422  {
423  const int j_E = d_indices[j_offset+e_j];
424  const int elt = j_E/elt_dofs;
425  j_elts[e_j] = elt;
426  j_B[e_j] = j_E%elt_dofs;
427  }
428  int min_e = GetMinElt(i_elts, i_nbElts, j_elts, j_nbElts);
429  if (e == min_e) // add the nnz only once
430  {
431  double val = 0.0;
432  for (int i = 0; i < i_nbElts; i++)
433  {
434  const int e_i = i_elts[i];
435  const int i_Bloc = i_B[i];
436  for (int j = 0; j < j_nbElts; j++)
437  {
438  const int e_j = j_elts[j];
439  const int j_Bloc = j_B[j];
440  if (e_i == e_j)
441  {
442  val += mat_ea(j_Bloc, i_Bloc, e_i);
443  }
444  }
445  }
446  const int nnz = GetAndIncrementNnzIndex(i_L, I);
447  J[nnz] = j_L;
448  Data[nnz] = val;
449  }
450  }
451  }
452  }
453  });
454  // We need to shift again the entries of I, we do it on CPU as it is very
455  // sequential.
456  auto h_I = mat.HostReadWriteI();
457  const int size = vd*all_dofs;
458  for (int i = 0; i < size; i++)
459  {
460  h_I[size-i] = h_I[size-(i+1)];
461  }
462  h_I[0] = 0;
463 }
464 
466  : ne(fes.GetNE()),
467  vdim(fes.GetVDim()),
468  byvdim(fes.GetOrdering() == Ordering::byVDIM),
469  ndof(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
470  ndofs(fes.GetNDofs())
471 {
472  height = vdim*ne*ndof;
473  width = vdim*ne*ndof;
474 }
475 
476 void L2ElementRestriction::Mult(const Vector &x, Vector &y) const
477 {
478  const int nd = ndof;
479  const int vd = vdim;
480  const bool t = byvdim;
481  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
482  auto d_y = Reshape(y.Write(), nd, vd, ne);
483  MFEM_FORALL(i, ndofs,
484  {
485  const int idx = i;
486  const int dof = idx % nd;
487  const int e = idx / nd;
488  for (int c = 0; c < vd; ++c)
489  {
490  d_y(dof, c, e) = d_x(t?c:idx, t?idx:c);
491  }
492  });
493 }
494 
496 {
497  const int nd = ndof;
498  const int vd = vdim;
499  const bool t = byvdim;
500  auto d_x = Reshape(x.Read(), nd, vd, ne);
501  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
502  MFEM_FORALL(i, ndofs,
503  {
504  const int idx = i;
505  const int dof = idx % nd;
506  const int e = idx / nd;
507  for (int c = 0; c < vd; ++c)
508  {
509  d_y(t?c:idx,t?idx:c) = d_x(dof, c, e);
510  }
511  });
512 }
513 
515 {
516  const int elem_dofs = ndof;
517  const int vd = vdim;
518  auto I = mat.WriteI();
519  const int isize = mat.Height() + 1;
520  const int interior_dofs = ne*elem_dofs*vd;
521  MFEM_FORALL(dof, isize,
522  {
523  I[dof] = dof<interior_dofs ? elem_dofs : 0;
524  });
525 }
526 
527 static MFEM_HOST_DEVICE int AddNnz(const int iE, int *I, const int dofs)
528 {
529  int val = AtomicAdd(I[iE],dofs);
530  return val;
531 }
532 
534  SparseMatrix &mat) const
535 {
536  const int elem_dofs = ndof;
537  const int vd = vdim;
538  auto I = mat.ReadWriteI();
539  auto J = mat.WriteJ();
540  auto Data = mat.WriteData();
541  auto mat_ea = Reshape(ea_data.Read(), elem_dofs, elem_dofs, ne);
542  MFEM_FORALL(iE, ne*elem_dofs*vd,
543  {
544  const int offset = AddNnz(iE,I,elem_dofs);
545  const int e = iE/elem_dofs;
546  const int i = iE%elem_dofs;
547  for (int j = 0; j < elem_dofs; j++)
548  {
549  J[offset+j] = e*elem_dofs+j;
550  Data[offset+j] = mat_ea(j,i,e);
551  }
552  });
553 }
554 
555 // Return the face degrees of freedom returned in Lexicographic order.
556 void GetFaceDofs(const int dim, const int face_id,
557  const int dof1d, Array<int> &faceMap)
558 {
559  switch (dim)
560  {
561  case 1:
562  switch (face_id)
563  {
564  case 0: // WEST
565  faceMap[0] = 0;
566  break;
567  case 1: // EAST
568  faceMap[0] = dof1d-1;
569  break;
570  }
571  break;
572  case 2:
573  switch (face_id)
574  {
575  case 0: // SOUTH
576  for (int i = 0; i < dof1d; ++i)
577  {
578  faceMap[i] = i;
579  }
580  break;
581  case 1: // EAST
582  for (int i = 0; i < dof1d; ++i)
583  {
584  faceMap[i] = dof1d-1 + i*dof1d;
585  }
586  break;
587  case 2: // NORTH
588  for (int i = 0; i < dof1d; ++i)
589  {
590  faceMap[i] = (dof1d-1)*dof1d + i;
591  }
592  break;
593  case 3: // WEST
594  for (int i = 0; i < dof1d; ++i)
595  {
596  faceMap[i] = i*dof1d;
597  }
598  break;
599  }
600  break;
601  case 3:
602  switch (face_id)
603  {
604  case 0: // BOTTOM
605  for (int i = 0; i < dof1d; ++i)
606  {
607  for (int j = 0; j < dof1d; ++j)
608  {
609  faceMap[i+j*dof1d] = i + j*dof1d;
610  }
611  }
612  break;
613  case 1: // SOUTH
614  for (int i = 0; i < dof1d; ++i)
615  {
616  for (int j = 0; j < dof1d; ++j)
617  {
618  faceMap[i+j*dof1d] = i + j*dof1d*dof1d;
619  }
620  }
621  break;
622  case 2: // EAST
623  for (int i = 0; i < dof1d; ++i)
624  {
625  for (int j = 0; j < dof1d; ++j)
626  {
627  faceMap[i+j*dof1d] = dof1d-1 + i*dof1d + j*dof1d*dof1d;
628  }
629  }
630  break;
631  case 3: // NORTH
632  for (int i = 0; i < dof1d; ++i)
633  {
634  for (int j = 0; j < dof1d; ++j)
635  {
636  faceMap[i+j*dof1d] = (dof1d-1)*dof1d + i + j*dof1d*dof1d;
637  }
638  }
639  break;
640  case 4: // WEST
641  for (int i = 0; i < dof1d; ++i)
642  {
643  for (int j = 0; j < dof1d; ++j)
644  {
645  faceMap[i+j*dof1d] = i*dof1d + j*dof1d*dof1d;
646  }
647  }
648  break;
649  case 5: // TOP
650  for (int i = 0; i < dof1d; ++i)
651  {
652  for (int j = 0; j < dof1d; ++j)
653  {
654  faceMap[i+j*dof1d] = (dof1d-1)*dof1d*dof1d + i + j*dof1d;
655  }
656  }
657  break;
658  }
659  break;
660  }
661 }
662 
664  const ElementDofOrdering e_ordering,
665  const FaceType type)
666  : fes(fes),
667  nf(fes.GetNFbyType(type)),
668  vdim(fes.GetVDim()),
669  byvdim(fes.GetOrdering() == Ordering::byVDIM),
670  ndofs(fes.GetNDofs()),
671  dof(nf > 0 ? fes.GetFaceElement(0)->GetDof() : 0),
672  nfdofs(nf*dof),
673  scatter_indices(nf*dof),
674  offsets(ndofs+1),
675  gather_indices(nf*dof)
676 {
677  if (nf==0) { return; }
678  // If fespace == H1
679  const FiniteElement *fe = fes.GetFE(0);
680  const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
681  MFEM_VERIFY(tfe != NULL &&
684  "Only Gauss-Lobatto and Bernstein basis are supported in "
685  "H1FaceRestriction.");
686  MFEM_VERIFY(fes.GetMesh()->Conforming(),
687  "Non-conforming meshes not yet supported with partial assembly.");
688  // Assuming all finite elements are using Gauss-Lobatto.
689  height = vdim*nf*dof;
690  width = fes.GetVSize();
691  const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
692  if (dof_reorder && nf > 0)
693  {
694  for (int f = 0; f < fes.GetNF(); ++f)
695  {
696  const FiniteElement *fe = fes.GetFaceElement(f);
697  const TensorBasisElement* el =
698  dynamic_cast<const TensorBasisElement*>(fe);
699  if (el) { continue; }
700  mfem_error("Finite element not suitable for lexicographic ordering");
701  }
702  const FiniteElement *fe = fes.GetFaceElement(0);
703  const TensorBasisElement* el =
704  dynamic_cast<const TensorBasisElement*>(fe);
705  const Array<int> &fe_dof_map = el->GetDofMap();
706  MFEM_VERIFY(fe_dof_map.Size() > 0, "invalid dof map");
707  }
708  const TensorBasisElement* el =
709  dynamic_cast<const TensorBasisElement*>(fe);
710  const int *dof_map = el->GetDofMap().GetData();
711  const Table& e2dTable = fes.GetElementToDofTable();
712  const int* elementMap = e2dTable.GetJ();
713  Array<int> faceMap(dof);
714  int e1, e2;
715  int inf1, inf2;
716  int face_id;
717  int orientation;
718  const int dof1d = fes.GetFE(0)->GetOrder()+1;
719  const int elem_dofs = fes.GetFE(0)->GetDof();
720  const int dim = fes.GetMesh()->SpaceDimension();
721  // Computation of scatter_indices
722  int f_ind = 0;
723  for (int f = 0; f < fes.GetNF(); ++f)
724  {
725  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
726  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
727  orientation = inf1 % 64;
728  face_id = inf1 / 64;
729  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
730  (type==FaceType::Boundary && e2<0 && inf2<0) )
731  {
732  // Assumes Gauss-Lobatto basis
733  if (dof_reorder)
734  {
735  if (orientation != 0)
736  {
737  mfem_error("FaceRestriction used on degenerated mesh.");
738  }
739  GetFaceDofs(dim, face_id, dof1d, faceMap); // Only for hex
740  }
741  else
742  {
743  mfem_error("FaceRestriction not yet implemented for this type of "
744  "element.");
745  // TODO Something with GetFaceDofs?
746  }
747  for (int d = 0; d < dof; ++d)
748  {
749  const int face_dof = faceMap[d];
750  const int did = (!dof_reorder)?face_dof:dof_map[face_dof];
751  const int gid = elementMap[e1*elem_dofs + did];
752  const int lid = dof*f_ind + d;
753  scatter_indices[lid] = gid;
754  }
755  f_ind++;
756  }
757  }
758  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
759  // Computation of gather_indices
760  for (int i = 0; i <= ndofs; ++i)
761  {
762  offsets[i] = 0;
763  }
764  f_ind = 0;
765  for (int f = 0; f < fes.GetNF(); ++f)
766  {
767  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
768  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
769  orientation = inf1 % 64;
770  face_id = inf1 / 64;
771  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
772  (type==FaceType::Boundary && e2<0 && inf2<0) )
773  {
774  GetFaceDofs(dim, face_id, dof1d, faceMap);
775  for (int d = 0; d < dof; ++d)
776  {
777  const int face_dof = faceMap[d];
778  const int did = (!dof_reorder)?face_dof:dof_map[face_dof];
779  const int gid = elementMap[e1*elem_dofs + did];
780  ++offsets[gid + 1];
781  }
782  f_ind++;
783  }
784  }
785  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
786  for (int i = 1; i <= ndofs; ++i)
787  {
788  offsets[i] += offsets[i - 1];
789  }
790  f_ind = 0;
791  for (int f = 0; f < fes.GetNF(); ++f)
792  {
793  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
794  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
795  orientation = inf1 % 64;
796  face_id = inf1 / 64;
797  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
798  (type==FaceType::Boundary && e2<0 && inf2<0) )
799  {
800  GetFaceDofs(dim, face_id, dof1d, faceMap);
801  for (int d = 0; d < dof; ++d)
802  {
803  const int face_dof = faceMap[d];
804  const int did = (!dof_reorder)?face_dof:dof_map[face_dof];
805  const int gid = elementMap[e1*elem_dofs + did];
806  const int lid = dof*f_ind + d;
807  gather_indices[offsets[gid]++] = lid;
808  }
809  f_ind++;
810  }
811  }
812  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
813  for (int i = ndofs; i > 0; --i)
814  {
815  offsets[i] = offsets[i - 1];
816  }
817  offsets[0] = 0;
818 }
819 
820 void H1FaceRestriction::Mult(const Vector& x, Vector& y) const
821 {
822  // Assumes all elements have the same number of dofs
823  const int nd = dof;
824  const int vd = vdim;
825  const bool t = byvdim;
826  auto d_indices = scatter_indices.Read();
827  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
828  auto d_y = Reshape(y.Write(), nd, vd, nf);
829  MFEM_FORALL(i, nfdofs,
830  {
831  const int idx = d_indices[i];
832  const int dof = i % nd;
833  const int face = i / nd;
834  for (int c = 0; c < vd; ++c)
835  {
836  d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
837  }
838  });
839 }
840 
842 {
843  // Assumes all elements have the same number of dofs
844  const int nd = dof;
845  const int vd = vdim;
846  const bool t = byvdim;
847  auto d_offsets = offsets.Read();
848  auto d_indices = gather_indices.Read();
849  auto d_x = Reshape(x.Read(), nd, vd, nf);
850  auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
851  MFEM_FORALL(i, ndofs,
852  {
853  const int offset = d_offsets[i];
854  const int nextOffset = d_offsets[i + 1];
855  for (int c = 0; c < vd; ++c)
856  {
857  double dofValue = 0;
858  for (int j = offset; j < nextOffset; ++j)
859  {
860  const int idx_j = d_indices[j];
861  dofValue += d_x(idx_j % nd, c, idx_j / nd);
862  }
863  d_y(t?c:i,t?i:c) += dofValue;
864  }
865  });
866 }
867 
868 static int ToLexOrdering2D(const int face_id, const int size1d, const int i)
869 {
870  if (face_id==2 || face_id==3)
871  {
872  return size1d-1-i;
873  }
874  else
875  {
876  return i;
877  }
878 }
879 
880 static int PermuteFace2D(const int face_id1, const int face_id2,
881  const int orientation,
882  const int size1d, const int index)
883 {
884  int new_index;
885  // Convert from lex ordering
886  if (face_id1==2 || face_id1==3)
887  {
888  new_index = size1d-1-index;
889  }
890  else
891  {
892  new_index = index;
893  }
894  // Permute based on face orientations
895  if (orientation==1)
896  {
897  new_index = size1d-1-new_index;
898  }
899  return ToLexOrdering2D(face_id2, size1d, new_index);
900 }
901 
902 static int ToLexOrdering3D(const int face_id, const int size1d, const int i,
903  const int j)
904 {
905  if (face_id==2 || face_id==1 || face_id==5)
906  {
907  return i + j*size1d;
908  }
909  else if (face_id==3 || face_id==4)
910  {
911  return (size1d-1-i) + j*size1d;
912  }
913  else // face_id==0
914  {
915  return i + (size1d-1-j)*size1d;
916  }
917 }
918 
919 static int PermuteFace3D(const int face_id1, const int face_id2,
920  const int orientation,
921  const int size1d, const int index)
922 {
923  int i=0, j=0, new_i=0, new_j=0;
924  i = index%size1d;
925  j = index/size1d;
926  // Convert from lex ordering
927  if (face_id1==3 || face_id1==4)
928  {
929  i = size1d-1-i;
930  }
931  else if (face_id1==0)
932  {
933  j = size1d-1-j;
934  }
935  // Permute based on face orientations
936  switch (orientation)
937  {
938  case 0:
939  new_i = i;
940  new_j = j;
941  break;
942  case 1:
943  new_i = j;
944  new_j = i;
945  break;
946  case 2:
947  new_i = j;
948  new_j = (size1d-1-i);
949  break;
950  case 3:
951  new_i = (size1d-1-i);
952  new_j = j;
953  break;
954  case 4:
955  new_i = (size1d-1-i);
956  new_j = (size1d-1-j);
957  break;
958  case 5:
959  new_i = (size1d-1-j);
960  new_j = (size1d-1-i);
961  break;
962  case 6:
963  new_i = (size1d-1-j);
964  new_j = i;
965  break;
966  case 7:
967  new_i = i;
968  new_j = (size1d-1-j);
969  break;
970  }
971  return ToLexOrdering3D(face_id2, size1d, new_i, new_j);
972 }
973 
974 // Permute dofs or quads on a face for e2 to match with the ordering of e1
975 int PermuteFaceL2(const int dim, const int face_id1,
976  const int face_id2, const int orientation,
977  const int size1d, const int index)
978 {
979  switch (dim)
980  {
981  case 1:
982  return 0;
983  case 2:
984  return PermuteFace2D(face_id1, face_id2, orientation, size1d, index);
985  case 3:
986  return PermuteFace3D(face_id1, face_id2, orientation, size1d, index);
987  default:
988  mfem_error("Unsupported dimension.");
989  return 0;
990  }
991 }
992 
994  const FaceType type,
995  const L2FaceValues m)
996  : fes(fes),
997  nf(fes.GetNFbyType(type)),
998  ne(fes.GetNE()),
999  vdim(fes.GetVDim()),
1000  byvdim(fes.GetOrdering() == Ordering::byVDIM),
1001  ndofs(fes.GetNDofs()),
1002  dof(nf > 0 ?
1003  fes.GetTraceElement(0, fes.GetMesh()->GetFaceBaseGeometry(0))->GetDof()
1004  : 0),
1005  elemDofs(fes.GetFE(0)->GetDof()),
1006  m(m),
1007  nfdofs(nf*dof),
1008  scatter_indices1(nf*dof),
1009  scatter_indices2(m==L2FaceValues::DoubleValued?nf*dof:0),
1010  offsets(ndofs+1),
1011  gather_indices((m==L2FaceValues::DoubleValued? 2 : 1)*nf*dof)
1012 {
1013 }
1014 
1016  const ElementDofOrdering e_ordering,
1017  const FaceType type,
1018  const L2FaceValues m)
1019  : L2FaceRestriction(fes, type, m)
1020 {
1021  // If fespace == L2
1022  const FiniteElement *fe = fes.GetFE(0);
1023  const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe);
1024  MFEM_VERIFY(tfe != NULL &&
1027  "Only Gauss-Lobatto and Bernstein basis are supported in "
1028  "L2FaceRestriction.");
1029  MFEM_VERIFY(fes.GetMesh()->Conforming(),
1030  "Non-conforming meshes not yet supported with partial assembly.");
1031  if (nf==0) { return; }
1033  width = fes.GetVSize();
1034  const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
1035  if (!dof_reorder)
1036  {
1037  mfem_error("Non-Tensor L2FaceRestriction not yet implemented.");
1038  }
1039  if (dof_reorder && nf > 0)
1040  {
1041  for (int f = 0; f < fes.GetNF(); ++f)
1042  {
1043  const FiniteElement *fe =
1045  const TensorBasisElement* el =
1046  dynamic_cast<const TensorBasisElement*>(fe);
1047  if (el) { continue; }
1048  mfem_error("Finite element not suitable for lexicographic ordering");
1049  }
1050  }
1051  const Table& e2dTable = fes.GetElementToDofTable();
1052  const int* elementMap = e2dTable.GetJ();
1053  Array<int> faceMap1(dof), faceMap2(dof);
1054  int e1, e2;
1055  int inf1, inf2;
1056  int face_id1 = -1, face_id2 = -1;
1057  int orientation = -1;
1058  const int dof1d = fes.GetFE(0)->GetOrder()+1;
1059  const int elem_dofs = fes.GetFE(0)->GetDof();
1060  const int dim = fes.GetMesh()->SpaceDimension();
1061  // Computation of scatter indices
1062  int f_ind=0;
1063  for (int f = 0; f < fes.GetNF(); ++f)
1064  {
1065  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
1066  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
1067  if (dof_reorder)
1068  {
1069  orientation = inf1 % 64;
1070  face_id1 = inf1 / 64;
1071  GetFaceDofs(dim, face_id1, dof1d, faceMap1); // Only for hex
1072  orientation = inf2 % 64;
1073  face_id2 = inf2 / 64;
1074  GetFaceDofs(dim, face_id2, dof1d, faceMap2); // Only for hex
1075  }
1076  else
1077  {
1078  mfem_error("FaceRestriction not yet implemented for this type of "
1079  "element.");
1080  // TODO Something with GetFaceDofs?
1081  }
1082  if ((type==FaceType::Interior && e2>=0) ||
1083  (type==FaceType::Boundary && e2<0))
1084  {
1085  for (int d = 0; d < dof; ++d)
1086  {
1087  const int face_dof = faceMap1[d];
1088  const int did = face_dof;
1089  const int gid = elementMap[e1*elem_dofs + did];
1090  const int lid = dof*f_ind + d;
1091  scatter_indices1[lid] = gid;
1092  }
1094  {
1095  for (int d = 0; d < dof; ++d)
1096  {
1097  if (type==FaceType::Interior && e2>=0) // interior face
1098  {
1099  const int pd = PermuteFaceL2(dim, face_id1, face_id2,
1100  orientation, dof1d, d);
1101  const int face_dof = faceMap2[pd];
1102  const int did = face_dof;
1103  const int gid = elementMap[e2*elem_dofs + did];
1104  const int lid = dof*f_ind + d;
1105  scatter_indices2[lid] = gid;
1106  }
1107  else if (type==FaceType::Boundary && e2<0) // true boundary face
1108  {
1109  const int lid = dof*f_ind + d;
1110  scatter_indices2[lid] = -1;
1111  }
1112  }
1113  }
1114  f_ind++;
1115  }
1116  }
1117  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
1118  // Computation of gather_indices
1119  for (int i = 0; i <= ndofs; ++i)
1120  {
1121  offsets[i] = 0;
1122  }
1123  f_ind = 0;
1124  for (int f = 0; f < fes.GetNF(); ++f)
1125  {
1126  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
1127  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
1128  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
1129  (type==FaceType::Boundary && e2<0 && inf2<0) )
1130  {
1131  orientation = inf1 % 64;
1132  face_id1 = inf1 / 64;
1133  GetFaceDofs(dim, face_id1, dof1d, faceMap1);
1134  orientation = inf2 % 64;
1135  face_id2 = inf2 / 64;
1136  GetFaceDofs(dim, face_id2, dof1d, faceMap2);
1137  for (int d = 0; d < dof; ++d)
1138  {
1139  const int did = faceMap1[d];
1140  const int gid = elementMap[e1*elem_dofs + did];
1141  ++offsets[gid + 1];
1142  }
1144  {
1145  for (int d = 0; d < dof; ++d)
1146  {
1147  if (type==FaceType::Interior && e2>=0) // interior face
1148  {
1149  const int pd = PermuteFaceL2(dim, face_id1, face_id2,
1150  orientation, dof1d, d);
1151  const int did = faceMap2[pd];
1152  const int gid = elementMap[e2*elem_dofs + did];
1153  ++offsets[gid + 1];
1154  }
1155  }
1156  }
1157  f_ind++;
1158  }
1159  }
1160  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
1161  for (int i = 1; i <= ndofs; ++i)
1162  {
1163  offsets[i] += offsets[i - 1];
1164  }
1165  f_ind = 0;
1166  for (int f = 0; f < fes.GetNF(); ++f)
1167  {
1168  fes.GetMesh()->GetFaceElements(f, &e1, &e2);
1169  fes.GetMesh()->GetFaceInfos(f, &inf1, &inf2);
1170  if ((type==FaceType::Interior && (e2>=0 || (e2<0 && inf2>=0))) ||
1171  (type==FaceType::Boundary && e2<0 && inf2<0) )
1172  {
1173  orientation = inf1 % 64;
1174  face_id1 = inf1 / 64;
1175  GetFaceDofs(dim, face_id1, dof1d, faceMap1);
1176  orientation = inf2 % 64;
1177  face_id2 = inf2 / 64;
1178  GetFaceDofs(dim, face_id2, dof1d, faceMap2);
1179  for (int d = 0; d < dof; ++d)
1180  {
1181  const int did = faceMap1[d];
1182  const int gid = elementMap[e1*elem_dofs + did];
1183  const int lid = dof*f_ind + d;
1184  // We don't shift lid to express that it's e1 of f
1185  gather_indices[offsets[gid]++] = lid;
1186  }
1188  {
1189  for (int d = 0; d < dof; ++d)
1190  {
1191  if (type==FaceType::Interior && e2>=0) // interior face
1192  {
1193  const int pd = PermuteFaceL2(dim, face_id1, face_id2,
1194  orientation, dof1d, d);
1195  const int did = faceMap2[pd];
1196  const int gid = elementMap[e2*elem_dofs + did];
1197  const int lid = dof*f_ind + d;
1198  // We shift lid to express that it's e2 of f
1199  gather_indices[offsets[gid]++] = nfdofs + lid;
1200  }
1201  }
1202  }
1203  f_ind++;
1204  }
1205  }
1206  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
1207  for (int i = ndofs; i > 0; --i)
1208  {
1209  offsets[i] = offsets[i - 1];
1210  }
1211  offsets[0] = 0;
1212 }
1213 
1214 void L2FaceRestriction::Mult(const Vector& x, Vector& y) const
1215 {
1216  // Assumes all elements have the same number of dofs
1217  const int nd = dof;
1218  const int vd = vdim;
1219  const bool t = byvdim;
1220 
1222  {
1223  auto d_indices1 = scatter_indices1.Read();
1224  auto d_indices2 = scatter_indices2.Read();
1225  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
1226  auto d_y = Reshape(y.Write(), nd, vd, 2, nf);
1227  MFEM_FORALL(i, nfdofs,
1228  {
1229  const int dof = i % nd;
1230  const int face = i / nd;
1231  const int idx1 = d_indices1[i];
1232  for (int c = 0; c < vd; ++c)
1233  {
1234  d_y(dof, c, 0, face) = d_x(t?c:idx1, t?idx1:c);
1235  }
1236  const int idx2 = d_indices2[i];
1237  for (int c = 0; c < vd; ++c)
1238  {
1239  d_y(dof, c, 1, face) = idx2==-1 ? 0.0 : d_x(t?c:idx2, t?idx2:c);
1240  }
1241  });
1242  }
1243  else
1244  {
1245  auto d_indices1 = scatter_indices1.Read();
1246  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
1247  auto d_y = Reshape(y.Write(), nd, vd, nf);
1248  MFEM_FORALL(i, nfdofs,
1249  {
1250  const int dof = i % nd;
1251  const int face = i / nd;
1252  const int idx1 = d_indices1[i];
1253  for (int c = 0; c < vd; ++c)
1254  {
1255  d_y(dof, c, face) = d_x(t?c:idx1, t?idx1:c);
1256  }
1257  });
1258  }
1259 }
1260 
1262 {
1263  // Assumes all elements have the same number of dofs
1264  const int nd = dof;
1265  const int vd = vdim;
1266  const bool t = byvdim;
1267  const int dofs = nfdofs;
1268  auto d_offsets = offsets.Read();
1269  auto d_indices = gather_indices.Read();
1270 
1272  {
1273  auto d_x = Reshape(x.Read(), nd, vd, 2, nf);
1274  auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
1275  MFEM_FORALL(i, ndofs,
1276  {
1277  const int offset = d_offsets[i];
1278  const int nextOffset = d_offsets[i + 1];
1279  for (int c = 0; c < vd; ++c)
1280  {
1281  double dofValue = 0;
1282  for (int j = offset; j < nextOffset; ++j)
1283  {
1284  int idx_j = d_indices[j];
1285  bool isE1 = idx_j < dofs;
1286  idx_j = isE1 ? idx_j : idx_j - dofs;
1287  dofValue += isE1 ?
1288  d_x(idx_j % nd, c, 0, idx_j / nd)
1289  :d_x(idx_j % nd, c, 1, idx_j / nd);
1290  }
1291  d_y(t?c:i,t?i:c) += dofValue;
1292  }
1293  });
1294  }
1295  else
1296  {
1297  auto d_x = Reshape(x.Read(), nd, vd, nf);
1298  auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
1299  MFEM_FORALL(i, ndofs,
1300  {
1301  const int offset = d_offsets[i];
1302  const int nextOffset = d_offsets[i + 1];
1303  for (int c = 0; c < vd; ++c)
1304  {
1305  double dofValue = 0;
1306  for (int j = offset; j < nextOffset; ++j)
1307  {
1308  int idx_j = d_indices[j];
1309  dofValue += d_x(idx_j % nd, c, idx_j / nd);
1310  }
1311  d_y(t?c:i,t?i:c) += dofValue;
1312  }
1313  });
1314  }
1315 }
1316 
1318  const bool keep_nbr_block) const
1319 {
1320  const int face_dofs = dof;
1321  auto d_indices1 = scatter_indices1.Read();
1322  auto d_indices2 = scatter_indices2.Read();
1323  auto I = mat.ReadWriteI();
1324  MFEM_FORALL(fdof, nf*face_dofs,
1325  {
1326  const int iE1 = d_indices1[fdof];
1327  const int iE2 = d_indices2[fdof];
1328  AddNnz(iE1,I,face_dofs);
1329  AddNnz(iE2,I,face_dofs);
1330  });
1331 }
1332 
1334  SparseMatrix &mat,
1335  const bool keep_nbr_block) const
1336 {
1337  const int face_dofs = dof;
1338  auto d_indices1 = scatter_indices1.Read();
1339  auto d_indices2 = scatter_indices2.Read();
1340  auto I = mat.ReadWriteI();
1341  auto mat_fea = Reshape(ea_data.Read(), face_dofs, face_dofs, 2, nf);
1342  auto J = mat.WriteJ();
1343  auto Data = mat.WriteData();
1344  MFEM_FORALL(fdof, nf*face_dofs,
1345  {
1346  const int f = fdof/face_dofs;
1347  const int iF = fdof%face_dofs;
1348  const int iE1 = d_indices1[f*face_dofs+iF];
1349  const int iE2 = d_indices2[f*face_dofs+iF];
1350  const int offset1 = AddNnz(iE1,I,face_dofs);
1351  const int offset2 = AddNnz(iE2,I,face_dofs);
1352  for (int jF = 0; jF < face_dofs; jF++)
1353  {
1354  const int jE1 = d_indices1[f*face_dofs+jF];
1355  const int jE2 = d_indices2[f*face_dofs+jF];
1356  J[offset2+jF] = jE1;
1357  J[offset1+jF] = jE2;
1358  Data[offset2+jF] = mat_fea(jF,iF,0,f);
1359  Data[offset1+jF] = mat_fea(jF,iF,1,f);
1360  }
1361  });
1362 }
1363 
1365  Vector &ea_data) const
1366 {
1367  const int face_dofs = dof;
1368  const int elem_dofs = elemDofs;
1369  const int NE = ne;
1371  {
1372  auto d_indices1 = scatter_indices1.Read();
1373  auto d_indices2 = scatter_indices2.Read();
1374  auto mat_fea = Reshape(fea_data.Read(), face_dofs, face_dofs, 2, nf);
1375  auto mat_ea = Reshape(ea_data.ReadWrite(), elem_dofs, elem_dofs, ne);
1376  MFEM_FORALL(f, nf,
1377  {
1378  const int e1 = d_indices1[f*face_dofs]/elem_dofs;
1379  const int e2 = d_indices2[f*face_dofs]/elem_dofs;
1380  for (int j = 0; j < face_dofs; j++)
1381  {
1382  const int jB1 = d_indices1[f*face_dofs+j]%elem_dofs;
1383  for (int i = 0; i < face_dofs; i++)
1384  {
1385  const int iB1 = d_indices1[f*face_dofs+i]%elem_dofs;
1386  AtomicAdd(mat_ea(iB1,jB1,e1), mat_fea(i,j,0,f));
1387  }
1388  }
1389  if (e2 < NE)
1390  {
1391  for (int j = 0; j < face_dofs; j++)
1392  {
1393  const int jB2 = d_indices2[f*face_dofs+j]%elem_dofs;
1394  for (int i = 0; i < face_dofs; i++)
1395  {
1396  const int iB2 = d_indices2[f*face_dofs+i]%elem_dofs;
1397  AtomicAdd(mat_ea(iB2,jB2,e2), mat_fea(i,j,1,f));
1398  }
1399  }
1400  }
1401  });
1402  }
1403  else
1404  {
1405  auto d_indices = scatter_indices1.Read();
1406  auto mat_fea = Reshape(fea_data.Read(), face_dofs, face_dofs, nf);
1407  auto mat_ea = Reshape(ea_data.ReadWrite(), elem_dofs, elem_dofs, ne);
1408  MFEM_FORALL(f, nf,
1409  {
1410  const int e = d_indices[f*face_dofs]/elem_dofs;
1411  for (int j = 0; j < face_dofs; j++)
1412  {
1413  const int jE = d_indices[f*face_dofs+j]%elem_dofs;
1414  for (int i = 0; i < face_dofs; i++)
1415  {
1416  const int iE = d_indices[f*face_dofs+i]%elem_dofs;
1417  AtomicAdd(mat_ea(iE,jE,e), mat_fea(i,j,f));
1418  }
1419  }
1420  });
1421  }
1422 }
1423 
1424 int ToLexOrdering(const int dim, const int face_id, const int size1d,
1425  const int index)
1426 {
1427  switch (dim)
1428  {
1429  case 1:
1430  return 0;
1431  case 2:
1432  return ToLexOrdering2D(face_id, size1d, index);
1433  case 3:
1434  return ToLexOrdering3D(face_id, size1d, index%size1d, index/size1d);
1435  default:
1436  mfem_error("Unsupported dimension.");
1437  return 0;
1438  }
1439 }
1440 
1441 } // namespace mfem
Abstract class for all finite elements.
Definition: fe.hpp:243
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition: array.hpp:320
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:540
void Mult(const Vector &x, Vector &y) const override
Extract the face degrees of freedom from x into y.
int * GetJ()
Definition: table.hpp:114
void MultTransposeUnsigned(const Vector &x, Vector &y) const
Compute MultTranspose without applying signs based on DOF orientations.
virtual void FillJAndData(const Vector &ea_data, SparseMatrix &mat, const bool keep_nbr_block=false) const
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:2759
bool Conforming() const
Definition: mesh.hpp:1366
Array< int > gather_indices
void AddMultTranspose(const Vector &x, Vector &y) const override
Add the face degrees of freedom x to the element degrees of freedom y.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition: array.hpp:304
MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
Definition: backends.hpp:62
T * GetData()
Returns the data.
Definition: array.hpp:108
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:327
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:438
void AddFaceMatricesToElementMatrices(Vector &fea_data, Vector &ea_data) const
This methods adds the DG face matrices to the element matrices.
void MultLeftInverse(const Vector &x, Vector &y) const
Geometry::Type GetFaceBaseGeometry(int i) const
Definition: mesh.hpp:968
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th face in the ...
Definition: fespace.cpp:2719
Memory< double > & GetMemoryData()
Definition: sparsemat.hpp:227
const L2FaceValues m
Memory< int > & GetMemoryI()
Definition: sparsemat.hpp:195
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
Array< int > gather_indices
L2ElementRestriction(const FiniteElementSpace &)
Memory< int > & GetMemoryJ()
Definition: sparsemat.hpp:211
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:136
int GetBasisType() const
Definition: fe.hpp:2146
FaceType
Definition: mesh.hpp:45
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
Bernstein polynomials.
Definition: fe.hpp:35
L2FaceRestriction(const FiniteElementSpace &, const FaceType, const L2FaceValues m=L2FaceValues::DoubleValued)
Data type sparse matrix.
Definition: sparsemat.hpp:41
int FillI(SparseMatrix &mat) const
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() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void FillI(SparseMatrix &mat) const
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:410
Operator that extracts Face degrees of freedom on L2 FiniteElementSpaces.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:434
void Mult(const Vector &x, Vector &y) const override
Extract the face degrees of freedom from x into y.
void BooleanMask(Vector &y) const
Fills the E-vector y with boolean values 0.0 and 1.0 such that each each entry of the L-vector is uni...
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
int * ReadWriteI(bool on_dev=true)
Definition: sparsemat.hpp:201
H1FaceRestriction(const FiniteElementSpace &fes, const ElementDofOrdering ordering, const FaceType type)
Constructor for a H1FaceRestriction.
void AddMultTranspose(const Vector &x, Vector &y) const override
Add the face degrees of freedom x to the element degrees of freedom y.
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:300
void MultUnsigned(const Vector &x, Vector &y) const
Compute Mult without applying signs based on DOF orientations.
int * WriteI(bool on_dev=true)
Definition: sparsemat.hpp:199
ElementRestriction(const FiniteElementSpace &, ElementDofOrdering)
Definition: restriction.cpp:21
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
Definition: fespace.hpp:573
int SpaceDimension() const
Definition: mesh.hpp:912
int ToLexOrdering(const int dim, const int face_id, const int size1d, const int index)
const FiniteElementSpace & fes
Definition: restriction.hpp:35
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1055
void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:1061
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:323
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
void GetFaceDofs(const int dim, const int face_id, const int dof1d, Array< int > &faceMap)
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition: fespace.hpp:28
int * WriteJ(bool on_dev=true)
Definition: sparsemat.hpp:215
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:2154
Array< int > scatter_indices1
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:442
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:65
int dim
Definition: ex24.cpp:53
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2388
int index(int i, int j, int nx, int ny)
Definition: life.cpp:237
Lexicographic ordering for tensor-product FiniteElements.
RefCoord t[3]
Vector data type.
Definition: vector.hpp:60
int PermuteFaceL2(const int dim, const int face_id1, const int face_id2, const int orientation, const int size1d, const int index)
const Table & GetElementToDofTable() const
Return a reference to the internal Table that stores the lists of scalar dofs, for each mesh element...
Definition: fespace.hpp:700
int * HostReadWriteI()
Definition: sparsemat.hpp:207
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:426
Array< int > scatter_indices
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
Array< int > scatter_indices2
virtual void FillI(SparseMatrix &mat, const bool keep_nbr_block=false) const
void FillSparseMatrix(const Vector &mat_ea, SparseMatrix &mat) const
Fill a Sparse Matrix with Element Matrices.
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
double * WriteData(bool on_dev=true)
Definition: sparsemat.hpp:231
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 ...