MFEM  v4.6.0
Finite element discretization library
restriction.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 #ifdef MFEM_USE_MPI
19 
20 #include "pfespace.hpp"
21 
22 #endif
23 
24 namespace mfem
25 {
26 
28  ElementDofOrdering e_ordering)
29  : fes(f),
30  ne(fes.GetNE()),
31  vdim(fes.GetVDim()),
32  byvdim(fes.GetOrdering() == Ordering::byVDIM),
33  ndofs(fes.GetNDofs()),
34  dof(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
35  nedofs(ne*dof),
36  offsets(ndofs+1),
37  indices(ne*dof),
38  gather_map(ne*dof)
39 {
40  // Assuming all finite elements are the same.
41  height = vdim*ne*dof;
42  width = fes.GetVSize();
43  const bool dof_reorder = (e_ordering == ElementDofOrdering::LEXICOGRAPHIC);
44  const int *dof_map = NULL;
45  if (dof_reorder && ne > 0)
46  {
47  for (int e = 0; e < ne; ++e)
48  {
49  const FiniteElement *fe = fes.GetFE(e);
50  const TensorBasisElement* el =
51  dynamic_cast<const TensorBasisElement*>(fe);
52  if (el) { continue; }
53  MFEM_ABORT("Finite element not suitable for lexicographic ordering");
54  }
55  const FiniteElement *fe = fes.GetFE(0);
56  const TensorBasisElement* el =
57  dynamic_cast<const TensorBasisElement*>(fe);
58  const Array<int> &fe_dof_map = el->GetDofMap();
59  MFEM_VERIFY(fe_dof_map.Size() > 0, "invalid dof map");
60  dof_map = fe_dof_map.GetData();
61  }
62  const Table& e2dTable = fes.GetElementToDofTable();
63  const int* element_map = e2dTable.GetJ();
64  // We will be keeping a count of how many local nodes point to its global dof
65  for (int i = 0; i <= ndofs; ++i)
66  {
67  offsets[i] = 0;
68  }
69  for (int e = 0; e < ne; ++e)
70  {
71  for (int d = 0; d < dof; ++d)
72  {
73  const int sgid = element_map[dof*e + d]; // signed
74  const int gid = (sgid >= 0) ? sgid : -1 - sgid;
75  ++offsets[gid + 1];
76  }
77  }
78  // Aggregate to find offsets for each global dof
79  for (int i = 1; i <= ndofs; ++i)
80  {
81  offsets[i] += offsets[i - 1];
82  }
83  // For each global dof, fill in all local nodes that point to it
84  for (int e = 0; e < ne; ++e)
85  {
86  for (int d = 0; d < dof; ++d)
87  {
88  const int sdid = dof_reorder ? dof_map[d] : 0; // signed
89  const int did = (!dof_reorder)?d:(sdid >= 0 ? sdid : -1-sdid);
90  const int sgid = element_map[dof*e + did]; // signed
91  const int gid = (sgid >= 0) ? sgid : -1-sgid;
92  const int lid = dof*e + d;
93  const bool plus = (sgid >= 0 && sdid >= 0) || (sgid < 0 && sdid < 0);
94  gather_map[lid] = plus ? gid : -1-gid;
95  indices[offsets[gid]++] = plus ? lid : -1-lid;
96  }
97  }
98  // We shifted the offsets vector by 1 by using it as a counter.
99  // Now we shift it back.
100  for (int i = ndofs; i > 0; --i)
101  {
102  offsets[i] = offsets[i - 1];
103  }
104  offsets[0] = 0;
105 }
106 
107 void ElementRestriction::Mult(const Vector& x, Vector& y) const
108 {
109  // Assumes all elements have the same number of dofs
110  const int nd = dof;
111  const int vd = vdim;
112  const bool t = byvdim;
113  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
114  auto d_y = Reshape(y.Write(), nd, vd, ne);
115  auto d_gather_map = gather_map.Read();
116  mfem::forall(dof*ne, [=] MFEM_HOST_DEVICE (int i)
117  {
118  const int gid = d_gather_map[i];
119  const bool plus = gid >= 0;
120  const int j = plus ? gid : -1-gid;
121  for (int c = 0; c < vd; ++c)
122  {
123  const double dof_value = d_x(t?c:j, t?j:c);
124  d_y(i % nd, c, i / nd) = plus ? dof_value : -dof_value;
125  }
126  });
127 }
128 
130 {
131  // Assumes all elements have the same number of dofs
132  const int nd = dof;
133  const int vd = vdim;
134  const bool t = byvdim;
135  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
136  auto d_y = Reshape(y.Write(), nd, vd, ne);
137  auto d_gather_map = gather_map.Read();
138 
139  mfem::forall(dof*ne, [=] MFEM_HOST_DEVICE (int i)
140  {
141  const int gid = d_gather_map[i];
142  const int j = gid >= 0 ? gid : -1-gid;
143  for (int c = 0; c < vd; ++c)
144  {
145  d_y(i % nd, c, i / nd) = d_x(t?c:j, t?j:c);
146  }
147  });
148 }
149 
150 template <bool ADD>
151 void ElementRestriction::TAddMultTranspose(const Vector& x, Vector& y) const
152 {
153  // Assumes all elements have the same number of dofs
154  const int nd = dof;
155  const int vd = vdim;
156  const bool t = byvdim;
157  auto d_offsets = offsets.Read();
158  auto d_indices = indices.Read();
159  auto d_x = Reshape(x.Read(), nd, vd, ne);
160  auto d_y = Reshape(ADD ? y.ReadWrite() : y.Write(), t?vd:ndofs, t?ndofs:vd);
161  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
162  {
163  const int offset = d_offsets[i];
164  const int next_offset = d_offsets[i + 1];
165  for (int c = 0; c < vd; ++c)
166  {
167  double dof_value = 0;
168  for (int j = offset; j < next_offset; ++j)
169  {
170  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
171  dof_value += ((d_indices[j] >= 0) ? d_x(idx_j % nd, c, idx_j / nd) :
172  -d_x(idx_j % nd, c, idx_j / nd));
173  }
174  if (ADD) { d_y(t?c:i,t?i:c) += dof_value; }
175  else { d_y(t?c:i,t?i:c) = dof_value; }
176  }
177  });
178 }
179 
181 {
182  constexpr bool ADD = false;
183  TAddMultTranspose<ADD>(x, y);
184 }
185 
187  const double a) const
188 {
189  MFEM_VERIFY(a == 1.0, "General coefficient case is not yet supported!");
190  constexpr bool ADD = true;
191  TAddMultTranspose<ADD>(x, y);
192 }
193 
195 {
196  // Assumes all elements have the same number of dofs
197  const int nd = dof;
198  const int vd = vdim;
199  const bool t = byvdim;
200  auto d_offsets = offsets.Read();
201  auto d_indices = indices.Read();
202  auto d_x = Reshape(x.Read(), nd, vd, ne);
203  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
204  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
205  {
206  const int offset = d_offsets[i];
207  const int next_offset = d_offsets[i + 1];
208  for (int c = 0; c < vd; ++c)
209  {
210  double dof_value = 0;
211  for (int j = offset; j < next_offset; ++j)
212  {
213  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
214  dof_value += d_x(idx_j % nd, c, idx_j / nd);
215  }
216  d_y(t?c:i,t?i:c) = dof_value;
217  }
218  });
219 }
220 
222 {
223  // Assumes all elements have the same number of dofs
224  const int nd = dof;
225  const int vd = vdim;
226  const bool t = byvdim;
227  auto d_offsets = offsets.Read();
228  auto d_indices = indices.Read();
229  auto d_x = Reshape(x.Read(), nd, vd, ne);
230  auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
231  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
232  {
233  const int next_offset = d_offsets[i + 1];
234  for (int c = 0; c < vd; ++c)
235  {
236  double dof_value = 0;
237  const int j = next_offset - 1;
238  const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
239  dof_value = (d_indices[j] >= 0) ? d_x(idx_j % nd, c, idx_j / nd) :
240  -d_x(idx_j % nd, c, idx_j / nd);
241  d_y(t?c:i,t?i:c) = dof_value;
242  }
243  });
244 }
245 
247 {
248  // Assumes all elements have the same number of dofs
249  const int nd = dof;
250  const int vd = vdim;
251  const bool t = byvdim;
252 
253  Array<char> processed(vd * ndofs);
254  processed = 0;
255 
256  auto d_offsets = offsets.HostRead();
257  auto d_indices = indices.HostRead();
258  auto d_x = Reshape(processed.HostReadWrite(), t?vd:ndofs, t?ndofs:vd);
259  auto d_y = Reshape(y.HostWrite(), nd, vd, ne);
260  for (int i = 0; i < ndofs; ++i)
261  {
262  const int offset = d_offsets[i];
263  const int next_offset = d_offsets[i+1];
264  for (int c = 0; c < vd; ++c)
265  {
266  for (int j = offset; j < next_offset; ++j)
267  {
268  const int idx_j = d_indices[j];
269  if (d_x(t?c:i,t?i:c))
270  {
271  d_y(idx_j % nd, c, idx_j / nd) = 0.0;
272  }
273  else
274  {
275  d_y(idx_j % nd, c, idx_j / nd) = 1.0;
276  d_x(t?c:i,t?i:c) = 1;
277  }
278  }
279  }
280  }
281 }
282 
284  SparseMatrix &mat) const
285 {
286  mat.GetMemoryI().New(mat.Height()+1, mat.GetMemoryI().GetMemoryType());
287  const int nnz = FillI(mat);
288  mat.GetMemoryJ().New(nnz, mat.GetMemoryJ().GetMemoryType());
289  mat.GetMemoryData().New(nnz, mat.GetMemoryData().GetMemoryType());
290  FillJAndData(mat_ea, mat);
291 }
292 
293 static MFEM_HOST_DEVICE int GetMinElt(const int *my_elts, const int nbElts,
294  const int *nbr_elts, const int nbrNbElts)
295 {
296  // Find the minimal element index found in both my_elts[] and nbr_elts[]
297  int min_el = INT_MAX;
298  for (int i = 0; i < nbElts; i++)
299  {
300  const int e_i = my_elts[i];
301  if (e_i >= min_el) { continue; }
302  for (int j = 0; j < nbrNbElts; j++)
303  {
304  if (e_i==nbr_elts[j])
305  {
306  min_el = e_i; // we already know e_i < min_el
307  break;
308  }
309  }
310  }
311  return min_el;
312 }
313 
314 /** Returns the index where a non-zero entry should be added and increment the
315  number of non-zeros for the row i_L. */
316 static MFEM_HOST_DEVICE int GetAndIncrementNnzIndex(const int i_L, int* I)
317 {
318  int ind = AtomicAdd(I[i_L],1);
319  return ind;
320 }
321 
323 {
324  static constexpr int Max = MaxNbNbr;
325  const int all_dofs = ndofs;
326  const int vd = vdim;
327  const int elt_dofs = dof;
328  auto I = mat.ReadWriteI();
329  auto d_offsets = offsets.Read();
330  auto d_indices = indices.Read();
331  auto d_gather_map = gather_map.Read();
332  mfem::forall(vd*all_dofs+1, [=] MFEM_HOST_DEVICE (int i_L)
333  {
334  I[i_L] = 0;
335  });
336  mfem::forall(ne*elt_dofs, [=] MFEM_HOST_DEVICE (int l_dof)
337  {
338  const int e = l_dof/elt_dofs;
339  const int i = l_dof%elt_dofs;
340 
341  int i_elts[Max];
342  const int i_gm = e*elt_dofs + i;
343  const int i_L = d_gather_map[i_gm];
344  const int i_offset = d_offsets[i_L];
345  const int i_next_offset = d_offsets[i_L+1];
346  const int i_nbElts = i_next_offset - i_offset;
347  MFEM_ASSERT_KERNEL(
348  i_nbElts <= Max,
349  "The connectivity of this mesh is beyond the max, increase the "
350  "MaxNbNbr variable to comply with your mesh.");
351  for (int e_i = 0; e_i < i_nbElts; ++e_i)
352  {
353  const int i_E = d_indices[i_offset+e_i];
354  i_elts[e_i] = i_E/elt_dofs;
355  }
356  for (int j = 0; j < elt_dofs; j++)
357  {
358  const int j_gm = e*elt_dofs + j;
359  const int j_L = d_gather_map[j_gm];
360  const int j_offset = d_offsets[j_L];
361  const int j_next_offset = d_offsets[j_L+1];
362  const int j_nbElts = j_next_offset - j_offset;
363  if (i_nbElts == 1 || j_nbElts == 1) // no assembly required
364  {
365  GetAndIncrementNnzIndex(i_L, I);
366  }
367  else // assembly required
368  {
369  int j_elts[Max];
370  for (int e_j = 0; e_j < j_nbElts; ++e_j)
371  {
372  const int j_E = d_indices[j_offset+e_j];
373  const int elt = j_E/elt_dofs;
374  j_elts[e_j] = elt;
375  }
376  int min_e = GetMinElt(i_elts, i_nbElts, j_elts, j_nbElts);
377  if (e == min_e) // add the nnz only once
378  {
379  GetAndIncrementNnzIndex(i_L, I);
380  }
381  }
382  }
383  });
384  // We need to sum the entries of I, we do it on CPU as it is very sequential.
385  auto h_I = mat.HostReadWriteI();
386  const int nTdofs = vd*all_dofs;
387  int sum = 0;
388  for (int i = 0; i < nTdofs; i++)
389  {
390  const int nnz = h_I[i];
391  h_I[i] = sum;
392  sum+=nnz;
393  }
394  h_I[nTdofs] = sum;
395  // We return the number of nnz
396  return h_I[nTdofs];
397 }
398 
400  SparseMatrix &mat) const
401 {
402  static constexpr int Max = MaxNbNbr;
403  const int all_dofs = ndofs;
404  const int vd = vdim;
405  const int elt_dofs = dof;
406  auto I = mat.ReadWriteI();
407  auto J = mat.WriteJ();
408  auto Data = mat.WriteData();
409  auto d_offsets = offsets.Read();
410  auto d_indices = indices.Read();
411  auto d_gather_map = gather_map.Read();
412  auto mat_ea = Reshape(ea_data.Read(), elt_dofs, elt_dofs, ne);
413  mfem::forall(ne*elt_dofs, [=] MFEM_HOST_DEVICE (int l_dof)
414  {
415  const int e = l_dof/elt_dofs;
416  const int i = l_dof%elt_dofs;
417 
418  int i_elts[Max];
419  int i_B[Max];
420  const int i_gm = e*elt_dofs + i;
421  const int i_L = d_gather_map[i_gm];
422  const int i_offset = d_offsets[i_L];
423  const int i_next_offset = d_offsets[i_L+1];
424  const int i_nbElts = i_next_offset - i_offset;
425  MFEM_ASSERT_KERNEL(
426  i_nbElts <= Max,
427  "The connectivity of this mesh is beyond the max, increase the "
428  "MaxNbNbr variable to comply with your mesh.");
429  for (int e_i = 0; e_i < i_nbElts; ++e_i)
430  {
431  const int i_E = d_indices[i_offset+e_i];
432  i_elts[e_i] = i_E/elt_dofs;
433  i_B[e_i] = i_E%elt_dofs;
434  }
435  for (int j = 0; j < elt_dofs; j++)
436  {
437  const int j_gm = e*elt_dofs + j;
438  const int j_L = d_gather_map[j_gm];
439  const int j_offset = d_offsets[j_L];
440  const int j_next_offset = d_offsets[j_L+1];
441  const int j_nbElts = j_next_offset - j_offset;
442  if (i_nbElts == 1 || j_nbElts == 1) // no assembly required
443  {
444  const int nnz = GetAndIncrementNnzIndex(i_L, I);
445  J[nnz] = j_L;
446  Data[nnz] = mat_ea(j,i,e);
447  }
448  else // assembly required
449  {
450  int j_elts[Max];
451  int j_B[Max];
452  for (int e_j = 0; e_j < j_nbElts; ++e_j)
453  {
454  const int j_E = d_indices[j_offset+e_j];
455  const int elt = j_E/elt_dofs;
456  j_elts[e_j] = elt;
457  j_B[e_j] = j_E%elt_dofs;
458  }
459  int min_e = GetMinElt(i_elts, i_nbElts, j_elts, j_nbElts);
460  if (e == min_e) // add the nnz only once
461  {
462  double val = 0.0;
463  for (int k = 0; k < i_nbElts; k++)
464  {
465  const int e_i = i_elts[k];
466  const int i_Bloc = i_B[k];
467  for (int l = 0; l < j_nbElts; l++)
468  {
469  const int e_j = j_elts[l];
470  const int j_Bloc = j_B[l];
471  if (e_i == e_j)
472  {
473  val += mat_ea(j_Bloc, i_Bloc, e_i);
474  }
475  }
476  }
477  const int nnz = GetAndIncrementNnzIndex(i_L, I);
478  J[nnz] = j_L;
479  Data[nnz] = val;
480  }
481  }
482  }
483  });
484  // We need to shift again the entries of I, we do it on CPU as it is very
485  // sequential.
486  auto h_I = mat.HostReadWriteI();
487  const int size = vd*all_dofs;
488  for (int i = 0; i < size; i++)
489  {
490  h_I[size-i] = h_I[size-(i+1)];
491  }
492  h_I[0] = 0;
493 }
494 
496  : ne(fes.GetNE()),
497  vdim(fes.GetVDim()),
498  byvdim(fes.GetOrdering() == Ordering::byVDIM),
499  ndof(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
500  ndofs(fes.GetNDofs())
501 {
502  height = vdim*ne*ndof;
503  width = vdim*ne*ndof;
504 }
505 
506 void L2ElementRestriction::Mult(const Vector &x, Vector &y) const
507 {
508  const int nd = ndof;
509  const int vd = vdim;
510  const bool t = byvdim;
511  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
512  auto d_y = Reshape(y.Write(), nd, vd, ne);
513  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
514  {
515  const int idx = i;
516  const int dof = idx % nd;
517  const int e = idx / nd;
518  for (int c = 0; c < vd; ++c)
519  {
520  d_y(dof, c, e) = d_x(t?c:idx, t?idx:c);
521  }
522  });
523 }
524 
525 template <bool ADD>
526 void L2ElementRestriction::TAddMultTranspose(const Vector &x, Vector &y) const
527 {
528  const int nd = ndof;
529  const int vd = vdim;
530  const bool t = byvdim;
531  auto d_x = Reshape(x.Read(), nd, vd, ne);
532  auto d_y = Reshape(ADD ? y.ReadWrite() : y.Write(), t?vd:ndofs, t?ndofs:vd);
533  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
534  {
535  const int idx = i;
536  const int dof = idx % nd;
537  const int e = idx / nd;
538  for (int c = 0; c < vd; ++c)
539  {
540  if (ADD) { d_y(t?c:idx,t?idx:c) += d_x(dof, c, e); }
541  else { d_y(t?c:idx,t?idx:c) = d_x(dof, c, e); }
542  }
543  });
544 }
545 
547 {
548  constexpr bool ADD = false;
549  TAddMultTranspose<ADD>(x, y);
550 }
551 
553  const double a) const
554 {
555  MFEM_VERIFY(a == 1.0, "General coefficient case is not yet supported!");
556  constexpr bool ADD = true;
557  TAddMultTranspose<ADD>(x, y);
558 }
559 
561 {
562  const int elem_dofs = ndof;
563  const int vd = vdim;
564  auto I = mat.WriteI();
565  const int isize = mat.Height() + 1;
566  const int interior_dofs = ne*elem_dofs*vd;
567  mfem::forall(isize, [=] MFEM_HOST_DEVICE (int dof)
568  {
569  I[dof] = dof<interior_dofs ? elem_dofs : 0;
570  });
571 }
572 
573 static MFEM_HOST_DEVICE int AddNnz(const int iE, int *I, const int dofs)
574 {
575  int val = AtomicAdd(I[iE],dofs);
576  return val;
577 }
578 
580  SparseMatrix &mat) const
581 {
582  const int elem_dofs = ndof;
583  const int vd = vdim;
584  auto I = mat.ReadWriteI();
585  auto J = mat.WriteJ();
586  auto Data = mat.WriteData();
587  auto mat_ea = Reshape(ea_data.Read(), elem_dofs, elem_dofs, ne);
588  mfem::forall(ne*elem_dofs*vd, [=] MFEM_HOST_DEVICE (int iE)
589  {
590  const int offset = AddNnz(iE,I,elem_dofs);
591  const int e = iE/elem_dofs;
592  const int i = iE%elem_dofs;
593  for (int j = 0; j < elem_dofs; j++)
594  {
595  J[offset+j] = e*elem_dofs+j;
596  Data[offset+j] = mat_ea(j,i,e);
597  }
598  });
599 }
600 
602  const FiniteElementSpace &fes,
603  const ElementDofOrdering f_ordering,
604  const FaceType type,
605  bool build)
606  : fes(fes),
607  nf(fes.GetNFbyType(type)),
608  vdim(fes.GetVDim()),
609  byvdim(fes.GetOrdering() == Ordering::byVDIM),
610  face_dofs(nf > 0 ? fes.GetFaceElement(0)->GetDof() : 0),
611  elem_dofs(fes.GetFE(0)->GetDof()),
612  nfdofs(nf*face_dofs),
613  ndofs(fes.GetNDofs()),
614  scatter_indices(nf*face_dofs),
615  gather_offsets(ndofs+1),
616  gather_indices(nf*face_dofs),
617  face_map(face_dofs)
618 {
620  width = fes.GetVSize();
621  if (nf==0) { return; }
622 
623  CheckFESpace(f_ordering);
624 
625  // Get the mapping from lexicographic DOF ordering to native ordering.
626  const TensorBasisElement* el =
627  dynamic_cast<const TensorBasisElement*>(fes.GetFE(0));
628  const Array<int> &dof_map_ = el->GetDofMap();
629  if (dof_map_.Size() > 0)
630  {
631  vol_dof_map.MakeRef(dof_map_);
632  }
633  else
634  {
635  // For certain types of elements dof_map_ is empty. In this case, that
636  // means the element is already ordered lexicographically, so the
637  // permutation is the identity.
639  for (int i = 0; i < elem_dofs; ++i) { vol_dof_map[i] = i; }
640  }
641 
642  if (!build) { return; }
643  ComputeScatterIndicesAndOffsets(f_ordering, type);
644  ComputeGatherIndices(f_ordering,type);
645 }
646 
648  const FiniteElementSpace &fes,
649  const ElementDofOrdering f_ordering,
650  const FaceType type)
651  : ConformingFaceRestriction(fes, f_ordering, type, true)
652 { }
653 
655 {
656  if (nf==0) { return; }
657  // Assumes all elements have the same number of dofs
658  const int nface_dofs = face_dofs;
659  const int vd = vdim;
660  const bool t = byvdim;
661  auto d_indices = scatter_indices.Read();
662  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
663  auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
664  mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
665  {
666  const int s_idx = d_indices[i];
667  const int sgn = (s_idx >= 0) ? 1 : -1;
668  const int idx = (s_idx >= 0) ? s_idx : -1 - s_idx;
669  const int dof = i % nface_dofs;
670  const int face = i / nface_dofs;
671  for (int c = 0; c < vd; ++c)
672  {
673  d_y(dof, c, face) = sgn*d_x(t?c:idx, t?idx:c);
674  }
675  });
676 }
677 
678 static void ConformingFaceRestriction_AddMultTranspose(
679  const int ndofs,
680  const int face_dofs,
681  const int nf,
682  const int vdim,
683  const bool by_vdim,
684  const Array<int> &gather_offsets,
685  const Array<int> &gather_indices,
686  const Vector &x,
687  Vector &y,
688  bool use_signs,
689  const double a)
690 {
691  MFEM_VERIFY(a == 1.0, "General coefficient case is not yet supported!");
692  if (nf==0) { return; }
693  // Assumes all elements have the same number of dofs
694  auto d_offsets = gather_offsets.Read();
695  auto d_indices = gather_indices.Read();
696  auto d_x = Reshape(x.Read(), face_dofs, vdim, nf);
697  auto d_y = Reshape(y.ReadWrite(), by_vdim?vdim:ndofs, by_vdim?ndofs:vdim);
698  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
699  {
700  const int offset = d_offsets[i];
701  const int next_offset = d_offsets[i + 1];
702  for (int c = 0; c < vdim; ++c)
703  {
704  double dof_value = 0;
705  for (int j = offset; j < next_offset; ++j)
706  {
707  const int s_idx_j = d_indices[j];
708  const double sgn = (s_idx_j >= 0 || !use_signs) ? 1.0 : -1.0;
709  const int idx_j = (s_idx_j >= 0) ? s_idx_j : -1 - s_idx_j;
710  dof_value += sgn*d_x(idx_j % face_dofs, c, idx_j / face_dofs);
711  }
712  d_y(by_vdim?c:i,by_vdim?i:c) += dof_value;
713  }
714  });
715 }
716 
718  const Vector& x, Vector& y, const double a) const
719 {
720  ConformingFaceRestriction_AddMultTranspose(
722  true, a);
723 }
724 
726  const Vector& x, Vector& y, const double a) const
727 {
728  ConformingFaceRestriction_AddMultTranspose(
730  false, a);
731 }
732 
734  f_ordering)
735 {
736 #ifdef MFEM_USE_MPI
737 
738  // If the underlying finite element space is parallel, ensure the face
739  // neighbor information is generated.
740  if (const ParFiniteElementSpace *pfes
741  = dynamic_cast<const ParFiniteElementSpace*>(&fes))
742  {
743  pfes->GetParMesh()->ExchangeFaceNbrData();
744  }
745 
746 #endif
747 
748 #ifdef MFEM_DEBUG
749  const FiniteElement *fe0 = fes.GetFE(0);
750  const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe0);
751  MFEM_VERIFY(tfe != NULL &&
754  "Only Gauss-Lobatto and Bernstein basis are supported in "
755  "ConformingFaceRestriction.");
756 
757  // Assuming all finite elements are using Gauss-Lobatto.
758  const bool dof_reorder = (f_ordering == ElementDofOrdering::LEXICOGRAPHIC);
759  if (dof_reorder && nf > 0)
760  {
761  for (int f = 0; f < fes.GetNF(); ++f)
762  {
763  const FiniteElement *fe = fes.GetFaceElement(f);
764  const TensorBasisElement* el =
765  dynamic_cast<const TensorBasisElement*>(fe);
766  if (el) { continue; }
767  MFEM_ABORT("Finite element not suitable for lexicographic ordering");
768  }
769  }
770 #endif
771 }
772 
773 void ConformingFaceRestriction::ComputeScatterIndicesAndOffsets(
774  const ElementDofOrdering f_ordering,
775  const FaceType type)
776 {
777  Mesh &mesh = *fes.GetMesh();
778 
779  // Initialization of the offsets
780  for (int i = 0; i <= ndofs; ++i)
781  {
782  gather_offsets[i] = 0;
783  }
784 
785  // Computation of scatter indices and offsets
786  int f_ind = 0;
787  for (int f = 0; f < fes.GetNF(); ++f)
788  {
789  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
790  if ( face.IsNonconformingCoarse() )
791  {
792  // We skip nonconforming coarse faces as they are treated
793  // by the corresponding nonconforming fine faces.
794  continue;
795  }
796  else if ( face.IsOfFaceType(type) )
797  {
798  SetFaceDofsScatterIndices(face, f_ind, f_ordering);
799  f_ind++;
800  }
801  }
802  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
803 
804  // Summation of the offsets
805  for (int i = 1; i <= ndofs; ++i)
806  {
807  gather_offsets[i] += gather_offsets[i - 1];
808  }
809 }
810 
811 void ConformingFaceRestriction::ComputeGatherIndices(
812  const ElementDofOrdering f_ordering,
813  const FaceType type)
814 {
815  Mesh &mesh = *fes.GetMesh();
816 
817  // Computation of gather_indices
818  int f_ind = 0;
819  for (int f = 0; f < fes.GetNF(); ++f)
820  {
821  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
822  if ( face.IsNonconformingCoarse() )
823  {
824  // We skip nonconforming coarse faces as they are treated
825  // by the corresponding nonconforming fine faces.
826  continue;
827  }
828  else if ( face.IsOfFaceType(type) )
829  {
830  SetFaceDofsGatherIndices(face, f_ind, f_ordering);
831  f_ind++;
832  }
833  }
834  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
835 
836  // Reset offsets to their initial value
837  for (int i = ndofs; i > 0; --i)
838  {
839  gather_offsets[i] = gather_offsets[i - 1];
840  }
841  gather_offsets[0] = 0;
842 }
843 
844 static inline int absdof(int i) { return i < 0 ? -1-i : i; }
845 
847  const Mesh::FaceInformation &face,
848  const int face_index,
849  const ElementDofOrdering f_ordering)
850 {
851  MFEM_ASSERT(!(face.IsNonconformingCoarse()),
852  "This method should not be used on nonconforming coarse faces.");
853  MFEM_ASSERT(face.element[0].orientation==0,
854  "FaceRestriction used on degenerated mesh.");
855  MFEM_CONTRACT_VAR(f_ordering); // not supported yet
856 
858 
859  const Table& e2dTable = fes.GetElementToDofTable();
860  const int* elem_map = e2dTable.GetJ();
861  const int elem_index = face.element[0].index;
862 
863  for (int face_dof = 0; face_dof < face_dofs; ++face_dof)
864  {
865  const int lex_volume_dof = face_map[face_dof];
866  const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof]; // signed
867  const int volume_dof = absdof(s_volume_dof);
868  const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
869  const int global_dof = absdof(s_global_dof);
870  const int restriction_dof = face_dofs*face_index + face_dof;
871  scatter_indices[restriction_dof] = s_global_dof;
872  ++gather_offsets[global_dof + 1];
873  }
874 }
875 
877  const Mesh::FaceInformation &face,
878  const int face_index,
879  const ElementDofOrdering f_ordering)
880 {
881  MFEM_ASSERT(!(face.IsNonconformingCoarse()),
882  "This method should not be used on nonconforming coarse faces.");
883  MFEM_CONTRACT_VAR(f_ordering); // not supported yet
884 
886 
887  const Table& e2dTable = fes.GetElementToDofTable();
888  const int* elem_map = e2dTable.GetJ();
889  const int elem_index = face.element[0].index;
890 
891  for (int face_dof = 0; face_dof < face_dofs; ++face_dof)
892  {
893  const int lex_volume_dof = face_map[face_dof];
894  const int s_volume_dof = AsConst(vol_dof_map)[lex_volume_dof];
895  const int volume_dof = absdof(s_volume_dof);
896  const int s_global_dof = elem_map[elem_index*elem_dofs + volume_dof];
897  const int sgn = (s_global_dof >= 0) ? 1 : -1;
898  const int global_dof = absdof(s_global_dof);
899  const int restriction_dof = face_dofs*face_index + face_dof;
900  const int s_restriction_dof = (sgn >= 0) ? restriction_dof : -1 -
901  restriction_dof;
902  gather_indices[gather_offsets[global_dof]++] = s_restriction_dof;
903  }
904 }
905 
906 static int ToLexOrdering2D(const int face_id, const int size1d, const int i)
907 {
908  if (face_id==2 || face_id==3)
909  {
910  return size1d-1-i;
911  }
912  else
913  {
914  return i;
915  }
916 }
917 
918 static int PermuteFace2D(const int face_id1, const int face_id2,
919  const int orientation,
920  const int size1d, const int index)
921 {
922  int new_index;
923  // Convert from lex ordering
924  if (face_id1==2 || face_id1==3)
925  {
926  new_index = size1d-1-index;
927  }
928  else
929  {
930  new_index = index;
931  }
932  // Permute based on face orientations
933  if (orientation==1)
934  {
935  new_index = size1d-1-new_index;
936  }
937  return ToLexOrdering2D(face_id2, size1d, new_index);
938 }
939 
940 static int ToLexOrdering3D(const int face_id, const int size1d, const int i,
941  const int j)
942 {
943  if (face_id==2 || face_id==1 || face_id==5)
944  {
945  return i + j*size1d;
946  }
947  else if (face_id==3 || face_id==4)
948  {
949  return (size1d-1-i) + j*size1d;
950  }
951  else // face_id==0
952  {
953  return i + (size1d-1-j)*size1d;
954  }
955 }
956 
957 static int PermuteFace3D(const int face_id1, const int face_id2,
958  const int orientation,
959  const int size1d, const int index)
960 {
961  int i=0, j=0, new_i=0, new_j=0;
962  i = index%size1d;
963  j = index/size1d;
964  // Convert from lex ordering
965  if (face_id1==3 || face_id1==4)
966  {
967  i = size1d-1-i;
968  }
969  else if (face_id1==0)
970  {
971  j = size1d-1-j;
972  }
973  // Permute based on face orientations
974  switch (orientation)
975  {
976  case 0:
977  new_i = i;
978  new_j = j;
979  break;
980  case 1:
981  new_i = j;
982  new_j = i;
983  break;
984  case 2:
985  new_i = j;
986  new_j = (size1d-1-i);
987  break;
988  case 3:
989  new_i = (size1d-1-i);
990  new_j = j;
991  break;
992  case 4:
993  new_i = (size1d-1-i);
994  new_j = (size1d-1-j);
995  break;
996  case 5:
997  new_i = (size1d-1-j);
998  new_j = (size1d-1-i);
999  break;
1000  case 6:
1001  new_i = (size1d-1-j);
1002  new_j = i;
1003  break;
1004  case 7:
1005  new_i = i;
1006  new_j = (size1d-1-j);
1007  break;
1008  }
1009  return ToLexOrdering3D(face_id2, size1d, new_i, new_j);
1010 }
1011 
1012 // Permute dofs or quads on a face for e2 to match with the ordering of e1
1013 int PermuteFaceL2(const int dim, const int face_id1,
1014  const int face_id2, const int orientation,
1015  const int size1d, const int index)
1016 {
1017  switch (dim)
1018  {
1019  case 1:
1020  return 0;
1021  case 2:
1022  return PermuteFace2D(face_id1, face_id2, orientation, size1d, index);
1023  case 3:
1024  return PermuteFace3D(face_id1, face_id2, orientation, size1d, index);
1025  default:
1026  MFEM_ABORT("Unsupported dimension.");
1027  return 0;
1028  }
1029 }
1030 
1032  const ElementDofOrdering f_ordering,
1033  const FaceType type,
1034  const L2FaceValues m,
1035  bool build)
1036  : fes(fes),
1037  nf(fes.GetNFbyType(type)),
1038  ne(fes.GetNE()),
1039  vdim(fes.GetVDim()),
1040  byvdim(fes.GetOrdering() == Ordering::byVDIM),
1041  face_dofs(nf > 0 ?
1042  fes.GetTraceElement(0, fes.GetMesh()->GetFaceGeometry(0))->GetDof()
1043  : 0),
1044  elem_dofs(fes.GetFE(0)->GetDof()),
1045  nfdofs(nf*face_dofs),
1046  ndofs(fes.GetNDofs()),
1047  type(type),
1048  m(m),
1049  scatter_indices1(nf*face_dofs),
1050  scatter_indices2(m==L2FaceValues::DoubleValued?nf*face_dofs:0),
1051  gather_offsets(ndofs+1),
1052  gather_indices((m==L2FaceValues::DoubleValued? 2 : 1)*nf*face_dofs),
1053  face_map(face_dofs)
1054 {
1056  width = fes.GetVSize();
1057  if (!build) { return; }
1058 
1059  CheckFESpace(f_ordering);
1060 
1061  ComputeScatterIndicesAndOffsets(f_ordering,type);
1062 
1063  ComputeGatherIndices(f_ordering, type);
1064 }
1065 
1067  const ElementDofOrdering f_ordering,
1068  const FaceType type,
1069  const L2FaceValues m)
1070  : L2FaceRestriction(fes, f_ordering, type, m, true)
1071 { }
1072 
1074  Vector& y) const
1075 {
1076  if (nf == 0) { return; }
1077  MFEM_ASSERT(
1079  "This method should be called when m == L2FaceValues::SingleValued.");
1080  // Assumes all elements have the same number of dofs
1081  const int nface_dofs = face_dofs;
1082  const int vd = vdim;
1083  const bool t = byvdim;
1084  auto d_indices1 = scatter_indices1.Read();
1085  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
1086  auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
1087  mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
1088  {
1089  const int dof = i % nface_dofs;
1090  const int face = i / nface_dofs;
1091  const int idx1 = d_indices1[i];
1092  for (int c = 0; c < vd; ++c)
1093  {
1094  d_y(dof, c, face) = d_x(t?c:idx1, t?idx1:c);
1095  }
1096  });
1097 }
1098 
1100  Vector& y) const
1101 {
1102  MFEM_ASSERT(
1104  "This method should be called when m == L2FaceValues::DoubleValued.");
1105  // Assumes all elements have the same number of dofs
1106  const int nface_dofs = face_dofs;
1107  const int vd = vdim;
1108  const bool t = byvdim;
1109  auto d_indices1 = scatter_indices1.Read();
1110  auto d_indices2 = scatter_indices2.Read();
1111  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
1112  auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
1113  mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
1114  {
1115  const int dof = i % nface_dofs;
1116  const int face = i / nface_dofs;
1117  const int idx1 = d_indices1[i];
1118  for (int c = 0; c < vd; ++c)
1119  {
1120  d_y(dof, c, 0, face) = d_x(t?c:idx1, t?idx1:c);
1121  }
1122  const int idx2 = d_indices2[i];
1123  for (int c = 0; c < vd; ++c)
1124  {
1125  d_y(dof, c, 1, face) = idx2==-1 ? 0.0 : d_x(t?c:idx2, t?idx2:c);
1126  }
1127  });
1128 }
1129 
1130 void L2FaceRestriction::Mult(const Vector& x, Vector& y) const
1131 {
1132  if (nf==0) { return; }
1134  {
1136  }
1137  else
1138  {
1140  }
1141 }
1142 
1144  const Vector& x, Vector& y) const
1145 {
1146  // Assumes all elements have the same number of dofs
1147  const int nface_dofs = face_dofs;
1148  const int vd = vdim;
1149  const bool t = byvdim;
1150  auto d_offsets = gather_offsets.Read();
1151  auto d_indices = gather_indices.Read();
1152  auto d_x = Reshape(x.Read(), nface_dofs, vd, nf);
1153  auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
1154  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
1155  {
1156  const int offset = d_offsets[i];
1157  const int next_offset = d_offsets[i + 1];
1158  for (int c = 0; c < vd; ++c)
1159  {
1160  double dof_value = 0;
1161  for (int j = offset; j < next_offset; ++j)
1162  {
1163  int idx_j = d_indices[j];
1164  dof_value += d_x(idx_j % nface_dofs, c, idx_j / nface_dofs);
1165  }
1166  d_y(t?c:i,t?i:c) += dof_value;
1167  }
1168  });
1169 }
1170 
1172  const Vector& x, Vector& y) const
1173 {
1174  // Assumes all elements have the same number of dofs
1175  const int nface_dofs = face_dofs;
1176  const int vd = vdim;
1177  const bool t = byvdim;
1178  const int dofs = nfdofs;
1179  auto d_offsets = gather_offsets.Read();
1180  auto d_indices = gather_indices.Read();
1181  auto d_x = Reshape(x.Read(), nface_dofs, vd, 2, nf);
1182  auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
1183  mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i)
1184  {
1185  const int offset = d_offsets[i];
1186  const int next_offset = d_offsets[i + 1];
1187  for (int c = 0; c < vd; ++c)
1188  {
1189  double dof_value = 0;
1190  for (int j = offset; j < next_offset; ++j)
1191  {
1192  int idx_j = d_indices[j];
1193  bool isE1 = idx_j < dofs;
1194  idx_j = isE1 ? idx_j : idx_j - dofs;
1195  dof_value += isE1 ?
1196  d_x(idx_j % nface_dofs, c, 0, idx_j / nface_dofs)
1197  :d_x(idx_j % nface_dofs, c, 1, idx_j / nface_dofs);
1198  }
1199  d_y(t?c:i,t?i:c) += dof_value;
1200  }
1201  });
1202 }
1203 
1205  const double a) const
1206 {
1207  MFEM_VERIFY(a == 1.0, "General coefficient case is not yet supported!");
1208  if (nf==0) { return; }
1210  {
1212  }
1213  else
1214  {
1216  }
1217 }
1218 
1220  const bool keep_nbr_block) const
1221 {
1222  const int nface_dofs = face_dofs;
1223  auto d_indices1 = scatter_indices1.Read();
1224  auto d_indices2 = scatter_indices2.Read();
1225  auto I = mat.ReadWriteI();
1226  mfem::forall(nf*nface_dofs, [=] MFEM_HOST_DEVICE (int fdof)
1227  {
1228  const int iE1 = d_indices1[fdof];
1229  const int iE2 = d_indices2[fdof];
1230  AddNnz(iE1,I,nface_dofs);
1231  AddNnz(iE2,I,nface_dofs);
1232  });
1233 }
1234 
1236  SparseMatrix &mat,
1237  const bool keep_nbr_block) const
1238 {
1239  const int nface_dofs = face_dofs;
1240  auto d_indices1 = scatter_indices1.Read();
1241  auto d_indices2 = scatter_indices2.Read();
1242  auto I = mat.ReadWriteI();
1243  auto mat_fea = Reshape(fea_data.Read(), nface_dofs, nface_dofs, 2, nf);
1244  auto J = mat.WriteJ();
1245  auto Data = mat.WriteData();
1246  mfem::forall(nf*nface_dofs, [=] MFEM_HOST_DEVICE (int fdof)
1247  {
1248  const int f = fdof/nface_dofs;
1249  const int iF = fdof%nface_dofs;
1250  const int iE1 = d_indices1[f*nface_dofs+iF];
1251  const int iE2 = d_indices2[f*nface_dofs+iF];
1252  const int offset1 = AddNnz(iE1,I,nface_dofs);
1253  const int offset2 = AddNnz(iE2,I,nface_dofs);
1254  for (int jF = 0; jF < nface_dofs; jF++)
1255  {
1256  const int jE1 = d_indices1[f*nface_dofs+jF];
1257  const int jE2 = d_indices2[f*nface_dofs+jF];
1258  J[offset2+jF] = jE1;
1259  J[offset1+jF] = jE2;
1260  Data[offset2+jF] = mat_fea(jF,iF,0,f);
1261  Data[offset1+jF] = mat_fea(jF,iF,1,f);
1262  }
1263  });
1264 }
1265 
1267  Vector &ea_data) const
1268 {
1269  const int nface_dofs = face_dofs;
1270  const int nelem_dofs = elem_dofs;
1271  const int NE = ne;
1273  {
1274  auto d_indices1 = scatter_indices1.Read();
1275  auto d_indices2 = scatter_indices2.Read();
1276  auto mat_fea = Reshape(fea_data.Read(), nface_dofs, nface_dofs, 2, nf);
1277  auto mat_ea = Reshape(ea_data.ReadWrite(), nelem_dofs, nelem_dofs, ne);
1278  mfem::forall(nf, [=] MFEM_HOST_DEVICE (int f)
1279  {
1280  const int e1 = d_indices1[f*nface_dofs]/nelem_dofs;
1281  const int e2 = d_indices2[f*nface_dofs]/nelem_dofs;
1282  for (int j = 0; j < nface_dofs; j++)
1283  {
1284  const int jB1 = d_indices1[f*nface_dofs+j]%nelem_dofs;
1285  for (int i = 0; i < nface_dofs; i++)
1286  {
1287  const int iB1 = d_indices1[f*nface_dofs+i]%nelem_dofs;
1288  AtomicAdd(mat_ea(iB1,jB1,e1), mat_fea(i,j,0,f));
1289  }
1290  }
1291  if (e2 < NE)
1292  {
1293  for (int j = 0; j < nface_dofs; j++)
1294  {
1295  const int jB2 = d_indices2[f*nface_dofs+j]%nelem_dofs;
1296  for (int i = 0; i < nface_dofs; i++)
1297  {
1298  const int iB2 = d_indices2[f*nface_dofs+i]%nelem_dofs;
1299  AtomicAdd(mat_ea(iB2,jB2,e2), mat_fea(i,j,1,f));
1300  }
1301  }
1302  }
1303  });
1304  }
1305  else
1306  {
1307  auto d_indices = scatter_indices1.Read();
1308  auto mat_fea = Reshape(fea_data.Read(), nface_dofs, nface_dofs, nf);
1309  auto mat_ea = Reshape(ea_data.ReadWrite(), nelem_dofs, nelem_dofs, ne);
1310  mfem::forall(nf, [=] MFEM_HOST_DEVICE (int f)
1311  {
1312  const int e = d_indices[f*nface_dofs]/nelem_dofs;
1313  for (int j = 0; j < nface_dofs; j++)
1314  {
1315  const int jE = d_indices[f*nface_dofs+j]%nelem_dofs;
1316  for (int i = 0; i < nface_dofs; i++)
1317  {
1318  const int iE = d_indices[f*nface_dofs+i]%nelem_dofs;
1319  AtomicAdd(mat_ea(iE,jE,e), mat_fea(i,j,f));
1320  }
1321  }
1322  });
1323  }
1324 }
1325 
1327 {
1328 #ifdef MFEM_USE_MPI
1329 
1330  // If the underlying finite element space is parallel, ensure the face
1331  // neighbor information is generated.
1332  if (const ParFiniteElementSpace *pfes
1333  = dynamic_cast<const ParFiniteElementSpace*>(&fes))
1334  {
1335  pfes->GetParMesh()->ExchangeFaceNbrData();
1336  }
1337 
1338 #endif
1339 
1340 #ifdef MFEM_DEBUG
1341  // If fespace == L2
1342  const FiniteElement *fe0 = fes.GetFE(0);
1343  const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe0);
1344  MFEM_VERIFY(tfe != NULL &&
1347  "Only Gauss-Lobatto and Bernstein basis are supported in "
1348  "L2FaceRestriction.");
1349  if (nf==0) { return; }
1350  const bool dof_reorder = (f_ordering == ElementDofOrdering::LEXICOGRAPHIC);
1351  if (!dof_reorder)
1352  {
1353  MFEM_ABORT("Non-Tensor L2FaceRestriction not yet implemented.");
1354  }
1355  if (dof_reorder && nf > 0)
1356  {
1357  for (int f = 0; f < fes.GetNF(); ++f)
1358  {
1359  const FiniteElement *fe =
1361  const TensorBasisElement* el =
1362  dynamic_cast<const TensorBasisElement*>(fe);
1363  if (el) { continue; }
1364  MFEM_ABORT("Finite element not suitable for lexicographic ordering");
1365  }
1366  }
1367 #endif
1368 }
1369 
1370 void L2FaceRestriction::ComputeScatterIndicesAndOffsets(
1371  const ElementDofOrdering f_ordering,
1372  const FaceType face_type)
1373 {
1374  Mesh &mesh = *fes.GetMesh();
1375  // Initialization of the offsets
1376  for (int i = 0; i <= ndofs; ++i)
1377  {
1378  gather_offsets[i] = 0;
1379  }
1380 
1381  // Computation of scatter indices and offsets
1382  int f_ind=0;
1383  for (int f = 0; f < fes.GetNF(); ++f)
1384  {
1385  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
1386  MFEM_ASSERT(!face.IsShared(),
1387  "Unexpected shared face in L2FaceRestriction.");
1388  if ( face.IsOfFaceType(face_type) )
1389  {
1390  SetFaceDofsScatterIndices1(face,f_ind);
1392  {
1393  if ( face_type==FaceType::Interior && face.IsInterior() )
1394  {
1396  }
1397  else if ( face_type==FaceType::Boundary && face.IsBoundary() )
1398  {
1399  SetBoundaryDofsScatterIndices2(face,f_ind);
1400  }
1401  }
1402  f_ind++;
1403  }
1404  }
1405  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
1406 
1407  // Summation of the offsets
1408  for (int i = 1; i <= ndofs; ++i)
1409  {
1410  gather_offsets[i] += gather_offsets[i - 1];
1411  }
1412 }
1413 
1414 void L2FaceRestriction::ComputeGatherIndices(
1415  const ElementDofOrdering f_ordering,
1416  const FaceType face_type)
1417 {
1418  Mesh &mesh = *fes.GetMesh();
1419  // Computation of gather_indices
1420  int f_ind = 0;
1421  for (int f = 0; f < fes.GetNF(); ++f)
1422  {
1423  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
1424  MFEM_ASSERT(!face.IsShared(),
1425  "Unexpected shared face in L2FaceRestriction.");
1426  if ( face.IsOfFaceType(face_type) )
1427  {
1428  SetFaceDofsGatherIndices1(face,f_ind);
1429  if ( m==L2FaceValues::DoubleValued &&
1430  face_type==FaceType::Interior &&
1431  face.IsLocal())
1432  {
1434  }
1435  f_ind++;
1436  }
1437  }
1438  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
1439 
1440  // Reset offsets to their correct value
1441  for (int i = ndofs; i > 0; --i)
1442  {
1443  gather_offsets[i] = gather_offsets[i - 1];
1444  }
1445  gather_offsets[0] = 0;
1446 }
1447 
1449  const Mesh::FaceInformation &face,
1450  const int face_index)
1451 {
1452  MFEM_ASSERT(!(face.IsNonconformingCoarse()),
1453  "This method should not be used on nonconforming coarse faces.");
1454  const Table& e2dTable = fes.GetElementToDofTable();
1455  const int* elem_map = e2dTable.GetJ();
1456  const int face_id1 = face.element[0].local_face_id;
1457  const int elem_index = face.element[0].index;
1458  fes.GetFE(0)->GetFaceMap(face_id1, face_map);
1459 
1460  for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
1461  {
1462  const int volume_dof_elem1 = face_map[face_dof_elem1];
1463  const int global_dof_elem1 = elem_map[elem_index*elem_dofs + volume_dof_elem1];
1464  const int restriction_dof_elem1 = face_dofs*face_index + face_dof_elem1;
1465  scatter_indices1[restriction_dof_elem1] = global_dof_elem1;
1466  ++gather_offsets[global_dof_elem1 + 1];
1467  }
1468 }
1469 
1471  const Mesh::FaceInformation &face,
1472  const int face_index)
1473 {
1474  MFEM_ASSERT(face.IsLocal(),
1475  "This method should only be used on local faces.");
1476  const Table& e2dTable = fes.GetElementToDofTable();
1477  const int* elem_map = e2dTable.GetJ();
1478  const int elem_index = face.element[1].index;
1479  const int face_id1 = face.element[0].local_face_id;
1480  const int face_id2 = face.element[1].local_face_id;
1481  const int orientation = face.element[1].orientation;
1482  const int dim = fes.GetMesh()->Dimension();
1483  const int dof1d = fes.GetFE(0)->GetOrder()+1;
1484  fes.GetFE(0)->GetFaceMap(face_id2, face_map);
1485 
1486  for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
1487  {
1488  const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
1489  orientation, dof1d,
1490  face_dof_elem1);
1491  const int volume_dof_elem2 = face_map[face_dof_elem2];
1492  const int global_dof_elem2 = elem_map[elem_index*elem_dofs + volume_dof_elem2];
1493  const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
1494  scatter_indices2[restriction_dof_elem2] = global_dof_elem2;
1495  ++gather_offsets[global_dof_elem2 + 1];
1496  }
1497 }
1498 
1500  const Mesh::FaceInformation &face,
1501  const int face_index)
1502 {
1503 #ifdef MFEM_USE_MPI
1504  MFEM_ASSERT(face.IsShared(),
1505  "This method should only be used on shared faces.");
1506  const int elem_index = face.element[1].index;
1507  const int face_id1 = face.element[0].local_face_id;
1508  const int face_id2 = face.element[1].local_face_id;
1509  const int orientation = face.element[1].orientation;
1510  const int dim = fes.GetMesh()->Dimension();
1511  const int dof1d = fes.GetFE(0)->GetOrder()+1;
1512  fes.GetFE(0)->GetFaceMap(face_id2, face_map);
1513  Array<int> face_nbr_dofs;
1514  const ParFiniteElementSpace &pfes =
1515  static_cast<const ParFiniteElementSpace&>(this->fes);
1516  pfes.GetFaceNbrElementVDofs(elem_index, face_nbr_dofs);
1517 
1518  for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
1519  {
1520  const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
1521  orientation, dof1d, face_dof_elem1);
1522  const int volume_dof_elem2 = face_map[face_dof_elem2];
1523  const int global_dof_elem2 = face_nbr_dofs[volume_dof_elem2];
1524  const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
1525  // Trick to differentiate dof location inter/shared
1526  scatter_indices2[restriction_dof_elem2] = ndofs+global_dof_elem2;
1527  }
1528 #endif
1529 }
1530 
1532  const Mesh::FaceInformation &face,
1533  const int face_index)
1534 {
1535  MFEM_ASSERT(face.IsBoundary(),
1536  "This method should only be used on boundary faces.");
1537 
1538  for (int d = 0; d < face_dofs; ++d)
1539  {
1540  const int restriction_dof_elem2 = face_dofs*face_index + d;
1541  scatter_indices2[restriction_dof_elem2] = -1;
1542  }
1543 }
1544 
1546  const Mesh::FaceInformation &face,
1547  const int face_index)
1548 {
1549  MFEM_ASSERT(!(face.IsNonconformingCoarse()),
1550  "This method should not be used on nonconforming coarse faces.");
1551  const Table& e2dTable = fes.GetElementToDofTable();
1552  const int* elem_map = e2dTable.GetJ();
1553  const int face_id1 = face.element[0].local_face_id;
1554  const int elem_index = face.element[0].index;
1555  fes.GetFE(0)->GetFaceMap(face_id1, face_map);
1556 
1557  for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
1558  {
1559  const int volume_dof_elem1 = face_map[face_dof_elem1];
1560  const int global_dof_elem1 = elem_map[elem_index*elem_dofs + volume_dof_elem1];
1561  const int restriction_dof_elem1 = face_dofs*face_index + face_dof_elem1;
1562  // We don't shift restriction_dof_elem1 to express that it's elem1 of the face
1563  gather_indices[gather_offsets[global_dof_elem1]++] = restriction_dof_elem1;
1564  }
1565 }
1566 
1568  const Mesh::FaceInformation &face,
1569  const int face_index)
1570 {
1571  MFEM_ASSERT(face.IsLocal(),
1572  "This method should only be used on local faces.");
1573  const Table& e2dTable = fes.GetElementToDofTable();
1574  const int* elem_map = e2dTable.GetJ();
1575  const int elem_index = face.element[1].index;
1576  const int face_id1 = face.element[0].local_face_id;
1577  const int face_id2 = face.element[1].local_face_id;
1578  const int orientation = face.element[1].orientation;
1579  const int dim = fes.GetMesh()->Dimension();
1580  const int dof1d = fes.GetFE(0)->GetOrder()+1;
1581  fes.GetFE(0)->GetFaceMap(face_id2, face_map);
1582 
1583  for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
1584  {
1585  const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
1586  orientation, dof1d,
1587  face_dof_elem1);
1588  const int volume_dof_elem2 = face_map[face_dof_elem2];
1589  const int global_dof_elem2 = elem_map[elem_index*elem_dofs + volume_dof_elem2];
1590  const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
1591  // We shift restriction_dof_elem2 to express that it's elem2 of the face
1592  gather_indices[gather_offsets[global_dof_elem2]++] = nfdofs +
1593  restriction_dof_elem2;
1594  }
1595 }
1596 
1598  ElementDofOrdering ordering,
1599  FaceType type)
1600  : fes(fes),
1601  ordering(ordering),
1602  interp_config(type==FaceType::Interior ? fes.GetNFbyType(type) : 0),
1603  nc_cpt(0)
1604 { }
1605 
1607  const Mesh::FaceInformation &face,
1608  int face_index)
1609 {
1610  interp_config[face_index] = InterpConfig();
1611 }
1612 
1614  const Mesh::FaceInformation &face,
1615  int face_index)
1616 {
1617  MFEM_ASSERT(!face.IsConforming(),
1618  "Registering face as nonconforming even though it is not.");
1619  const DenseMatrix* ptMat = face.point_matrix;
1620  // In the case of nonconforming slave shared face the master face is elem1.
1621  const int master_side =
1623  const int face_key = (master_side == 0 ? 1000 : 0) +
1624  face.element[0].local_face_id +
1625  6*face.element[1].local_face_id +
1626  36*face.element[1].orientation ;
1627  // Unfortunately we can't trust unicity of the ptMat to identify the transformation.
1628  Key key(ptMat, face_key);
1629  auto itr = interp_map.find(key);
1630  if ( itr == interp_map.end() )
1631  {
1632  const DenseMatrix* interpolator =
1633  GetCoarseToFineInterpolation(face,ptMat);
1634  interp_map[key] = {nc_cpt, interpolator};
1635  interp_config[face_index] = {master_side, nc_cpt};
1636  nc_cpt++;
1637  }
1638  else
1639  {
1640  interp_config[face_index] = {master_side, itr->second.first};
1641  }
1642 }
1643 
1644 const DenseMatrix* InterpolationManager::GetCoarseToFineInterpolation(
1645  const Mesh::FaceInformation &face,
1646  const DenseMatrix* ptMat)
1647 {
1649  "The following interpolation operator is only implemented for"
1650  "lexicographic ordering.");
1651  MFEM_VERIFY(!face.IsConforming(),
1652  "This method should not be called on conforming faces.")
1653  const int face_id1 = face.element[0].local_face_id;
1654  const int face_id2 = face.element[1].local_face_id;
1655 
1656  const bool is_ghost_slave =
1658  const int master_face_id = is_ghost_slave ? face_id1 : face_id2;
1659 
1660  // Computation of the interpolation matrix from master
1661  // (coarse) face to slave (fine) face.
1662  // Assumes all trace elements are the same.
1663  const FiniteElement *trace_fe =
1665  const int face_dofs = trace_fe->GetDof();
1666  const TensorBasisElement* el =
1667  dynamic_cast<const TensorBasisElement*>(trace_fe);
1668  const auto dof_map = el->GetDofMap();
1669  DenseMatrix* interpolator = new DenseMatrix(face_dofs,face_dofs);
1670  Vector shape(face_dofs);
1671 
1673  isotr.SetIdentityTransformation(trace_fe->GetGeomType());
1674  isotr.SetPointMat(*ptMat);
1675  DenseMatrix& trans_pt_mat = isotr.GetPointMat();
1676  // PointMatrix needs to be flipped in 2D
1677  if ( trace_fe->GetGeomType()==Geometry::SEGMENT && !is_ghost_slave )
1678  {
1679  std::swap(trans_pt_mat(0,0),trans_pt_mat(0,1));
1680  }
1681  DenseMatrix native_interpolator(face_dofs,face_dofs);
1682  trace_fe->GetLocalInterpolation(isotr, native_interpolator);
1683  const int dim = trace_fe->GetDim()+1;
1684  const int dof1d = trace_fe->GetOrder()+1;
1685  const int orientation = face.element[1].orientation;
1686  for (int i = 0; i < face_dofs; i++)
1687  {
1688  const int ni = (dof_map.Size()==0) ? i : dof_map[i];
1689  int li = ToLexOrdering(dim, master_face_id, dof1d, i);
1690  if ( !is_ghost_slave )
1691  {
1692  // master side is elem 2, so we permute to order dofs as elem 1.
1693  li = PermuteFaceL2(dim, face_id2, face_id1,
1694  orientation, dof1d, li);
1695  }
1696  for (int j = 0; j < face_dofs; j++)
1697  {
1698  int lj = ToLexOrdering(dim, master_face_id, dof1d, j);
1699  if ( !is_ghost_slave )
1700  {
1701  // master side is elem 2, so we permute to order dofs as elem 1.
1702  lj = PermuteFaceL2(dim, face_id2, face_id1,
1703  orientation, dof1d, lj);
1704  }
1705  const int nj = (dof_map.Size()==0) ? j : dof_map[j];
1706  (*interpolator)(li,lj) = native_interpolator(ni,nj);
1707  }
1708  }
1709  return interpolator;
1710 }
1711 
1713 {
1714  // Assumes all trace elements are the same.
1715  const FiniteElement *trace_fe =
1717  const int face_dofs = trace_fe->GetDof();
1718  const int nc_size = interp_map.size();
1719  MFEM_VERIFY(nc_cpt==nc_size, "Unexpected number of interpolators.");
1720  interpolators.SetSize(face_dofs*face_dofs*nc_size);
1721  auto d_interp = Reshape(interpolators.HostWrite(),face_dofs,face_dofs,nc_size);
1722  for (auto val : interp_map)
1723  {
1724  const int idx = val.second.first;
1725  const DenseMatrix &interpolator = *val.second.second;
1726  for (int i = 0; i < face_dofs; i++)
1727  {
1728  for (int j = 0; j < face_dofs; j++)
1729  {
1730  d_interp(i,j,idx) = interpolator(i,j);
1731  }
1732  }
1733  delete val.second.second;
1734  }
1735  interp_map.clear();
1736 }
1737 
1739 {
1740  // Count nonconforming faces
1741  int num_nc_faces = 0;
1742  for (int i = 0; i < interp_config.Size(); i++)
1743  {
1744  if ( interp_config[i].is_non_conforming )
1745  {
1746  num_nc_faces++;
1747  }
1748  }
1749  // Set nc_interp_config
1750  nc_interp_config.SetSize(num_nc_faces);
1751  int nc_index = 0;
1752  for (int i = 0; i < interp_config.Size(); i++)
1753  {
1754  auto & config = interp_config[i];
1755  if ( config.is_non_conforming )
1756  {
1757  nc_interp_config[nc_index] = NCInterpConfig(i, config);
1758  nc_index++;
1759  }
1760  }
1761 }
1762 
1764  const ElementDofOrdering f_ordering,
1765  const FaceType type,
1766  const L2FaceValues m,
1767  bool build)
1768  : L2FaceRestriction(fes, f_ordering, type, m, false),
1769  interpolations(fes, f_ordering, type)
1770 {
1771  if (!build) { return; }
1772  x_interp.UseDevice(true);
1773 
1774  CheckFESpace(f_ordering);
1775 
1776  ComputeScatterIndicesAndOffsets(f_ordering, type);
1777 
1778  ComputeGatherIndices(f_ordering, type);
1779 }
1780 
1782  const ElementDofOrdering f_ordering,
1783  const FaceType type,
1784  const L2FaceValues m)
1785  : NCL2FaceRestriction(fes, f_ordering, type, m, true)
1786 { }
1787 
1789  const Vector& x, Vector& y) const
1790 {
1793 }
1794 
1796  Vector& y) const
1797 {
1798  if (nf == 0) { return; }
1799  // Assumes all elements have the same number of dofs
1800  const int nface_dofs = face_dofs;
1801  const int vd = vdim;
1802  auto d_y = Reshape(y.ReadWrite(), nface_dofs, vd, 2, nf);
1803  auto &nc_interp_config = interpolations.GetNCFaceInterpConfig();
1804  const int num_nc_faces = nc_interp_config.Size();
1805  if ( num_nc_faces == 0 ) { return; }
1806  auto interp_config_ptr = nc_interp_config.Read();
1807  const int nc_size = interpolations.GetNumInterpolators();
1808  auto d_interp = Reshape(interpolations.GetInterpolators().Read(),
1809  nface_dofs, nface_dofs, nc_size);
1810  static constexpr int max_nd = 16*16;
1811  MFEM_VERIFY(nface_dofs<=max_nd, "Too many degrees of freedom.");
1812  mfem::forall_2D(num_nc_faces, nface_dofs, 1, [=] MFEM_HOST_DEVICE (int nc_face)
1813  {
1814  MFEM_SHARED double dof_values[max_nd];
1815  const NCInterpConfig conf = interp_config_ptr[nc_face];
1816  if ( conf.is_non_conforming )
1817  {
1818  const int master_side = conf.master_side;
1819  const int interp_index = conf.index;
1820  const int face = conf.face_index;
1821  for (int c = 0; c < vd; ++c)
1822  {
1823  MFEM_FOREACH_THREAD(dof,x,nface_dofs)
1824  {
1825  dof_values[dof] = d_y(dof, c, master_side, face);
1826  }
1827  MFEM_SYNC_THREAD;
1828  MFEM_FOREACH_THREAD(dof_out,x,nface_dofs)
1829  {
1830  double res = 0.0;
1831  for (int dof_in = 0; dof_in<nface_dofs; dof_in++)
1832  {
1833  res += d_interp(dof_out, dof_in, interp_index)*dof_values[dof_in];
1834  }
1835  d_y(dof_out, c, master_side, face) = res;
1836  }
1837  MFEM_SYNC_THREAD;
1838  }
1839  }
1840  });
1841 }
1842 
1843 void NCL2FaceRestriction::Mult(const Vector& x, Vector& y) const
1844 {
1845  if ( type==FaceType::Interior && m==L2FaceValues::DoubleValued )
1846  {
1847  DoubleValuedNonconformingMult(x, y);
1848  }
1849  else if ( type==FaceType::Boundary && m==L2FaceValues::DoubleValued )
1850  {
1851  DoubleValuedConformingMult(x, y);
1852  }
1853  else // Single valued (assumes no nonconforming master on elem1)
1854  {
1855  SingleValuedConformingMult(x, y);
1856  }
1857 }
1858 
1859 void NCL2FaceRestriction::SingleValuedNonconformingTransposeInterpolation(
1860  const Vector& x) const
1861 {
1862  MFEM_ASSERT(
1863  m == L2FaceValues::SingleValued,
1864  "This method should be called when m == L2FaceValues::SingleValued.");
1865  if (x_interp.Size()==0)
1866  {
1867  x_interp.SetSize(x.Size());
1868  }
1869  x_interp = x;
1870  SingleValuedNonconformingTransposeInterpolationInPlace(x_interp);
1871 }
1872 
1873 
1874 void NCL2FaceRestriction::SingleValuedNonconformingTransposeInterpolationInPlace(
1875  Vector& x) const
1876 {
1877  // Assumes all elements have the same number of dofs
1878  const int nface_dofs = face_dofs;
1879  const int vd = vdim;
1880  // Interpolation
1881  auto d_x = Reshape(x_interp.ReadWrite(), nface_dofs, vd, nf);
1882  auto &nc_interp_config = interpolations.GetNCFaceInterpConfig();
1883  const int num_nc_faces = nc_interp_config.Size();
1884  if ( num_nc_faces == 0 ) { return; }
1885  auto interp_config_ptr = nc_interp_config.Read();
1886  auto interpolators = interpolations.GetInterpolators().Read();
1887  const int nc_size = interpolations.GetNumInterpolators();
1888  auto d_interp = Reshape(interpolators, nface_dofs, nface_dofs, nc_size);
1889  static constexpr int max_nd = 16*16;
1890  MFEM_VERIFY(nface_dofs<=max_nd, "Too many degrees of freedom.");
1891  mfem::forall_2D(num_nc_faces, nface_dofs, 1, [=] MFEM_HOST_DEVICE (int nc_face)
1892  {
1893  MFEM_SHARED double dof_values[max_nd];
1894  const NCInterpConfig conf = interp_config_ptr[nc_face];
1895  const int master_side = conf.master_side;
1896  const int interp_index = conf.index;
1897  const int face = conf.face_index;
1898  if ( conf.is_non_conforming && master_side==0 )
1899  {
1900  // Interpolation from fine to coarse
1901  for (int c = 0; c < vd; ++c)
1902  {
1903  MFEM_FOREACH_THREAD(dof,x,nface_dofs)
1904  {
1905  dof_values[dof] = d_x(dof, c, face);
1906  }
1907  MFEM_SYNC_THREAD;
1908  MFEM_FOREACH_THREAD(dof_out,x,nface_dofs)
1909  {
1910  double res = 0.0;
1911  for (int dof_in = 0; dof_in<nface_dofs; dof_in++)
1912  {
1913  res += d_interp(dof_in, dof_out, interp_index)*dof_values[dof_in];
1914  }
1915  d_x(dof_out, c, face) = res;
1916  }
1917  MFEM_SYNC_THREAD;
1918  }
1919  }
1920  });
1921 }
1922 
1923 void NCL2FaceRestriction::DoubleValuedNonconformingTransposeInterpolation(
1924  const Vector& x) const
1925 {
1926  MFEM_ASSERT(
1927  m == L2FaceValues::DoubleValued,
1928  "This method should be called when m == L2FaceValues::DoubleValued.");
1929  if (x_interp.Size()==0)
1930  {
1931  x_interp.SetSize(x.Size());
1932  }
1933  x_interp = x;
1934  DoubleValuedNonconformingTransposeInterpolationInPlace(x_interp);
1935 }
1936 
1937 void NCL2FaceRestriction::DoubleValuedNonconformingTransposeInterpolationInPlace(
1938  Vector& x) const
1939 {
1940  // Assumes all elements have the same number of dofs
1941  const int nface_dofs = face_dofs;
1942  const int vd = vdim;
1943  // Interpolation
1944  auto d_x = Reshape(x.ReadWrite(), nface_dofs, vd, 2, nf);
1945  auto &nc_interp_config = interpolations.GetNCFaceInterpConfig();
1946  const int num_nc_faces = nc_interp_config.Size();
1947  if ( num_nc_faces == 0 ) { return; }
1948  auto interp_config_ptr = nc_interp_config.Read();
1949  auto interpolators = interpolations.GetInterpolators().Read();
1950  const int nc_size = interpolations.GetNumInterpolators();
1951  auto d_interp = Reshape(interpolators, nface_dofs, nface_dofs, nc_size);
1952  static constexpr int max_nd = 16*16;
1953  MFEM_VERIFY(nface_dofs<=max_nd, "Too many degrees of freedom.");
1954  mfem::forall_2D(num_nc_faces, nface_dofs, 1, [=] MFEM_HOST_DEVICE (int nc_face)
1955  {
1956  MFEM_SHARED double dof_values[max_nd];
1957  const NCInterpConfig conf = interp_config_ptr[nc_face];
1958  const int master_side = conf.master_side;
1959  const int interp_index = conf.index;
1960  const int face = conf.face_index;
1961  if ( conf.is_non_conforming )
1962  {
1963  // Interpolation from fine to coarse
1964  for (int c = 0; c < vd; ++c)
1965  {
1966  MFEM_FOREACH_THREAD(dof,x,nface_dofs)
1967  {
1968  dof_values[dof] = d_x(dof, c, master_side, face);
1969  }
1970  MFEM_SYNC_THREAD;
1971  MFEM_FOREACH_THREAD(dof_out,x,nface_dofs)
1972  {
1973  double res = 0.0;
1974  for (int dof_in = 0; dof_in<nface_dofs; dof_in++)
1975  {
1976  res += d_interp(dof_in, dof_out, interp_index)*dof_values[dof_in];
1977  }
1978  d_x(dof_out, c, master_side, face) = res;
1979  }
1980  MFEM_SYNC_THREAD;
1981  }
1982  }
1983  });
1984 }
1985 
1986 void NCL2FaceRestriction::AddMultTranspose(const Vector& x, Vector& y,
1987  const double a) const
1988 {
1989  MFEM_VERIFY(a == 1.0, "General coefficient case is not yet supported!");
1990  if (nf==0) { return; }
1991  if (type==FaceType::Interior)
1992  {
1993  if ( m==L2FaceValues::DoubleValued )
1994  {
1995  DoubleValuedNonconformingTransposeInterpolation(x);
1996  DoubleValuedConformingAddMultTranspose(x_interp, y);
1997  }
1998  else if ( m==L2FaceValues::SingleValued )
1999  {
2000  SingleValuedNonconformingTransposeInterpolation(x);
2001  SingleValuedConformingAddMultTranspose(x_interp, y);
2002  }
2003  }
2004  else
2005  {
2006  if ( m==L2FaceValues::DoubleValued )
2007  {
2008  DoubleValuedConformingAddMultTranspose(x, y);
2009  }
2010  else if ( m==L2FaceValues::SingleValued )
2011  {
2012  SingleValuedConformingAddMultTranspose(x, y);
2013  }
2014  }
2015 }
2016 
2017 void NCL2FaceRestriction::AddMultTransposeInPlace(Vector& x, Vector& y) const
2018 {
2019  if (nf==0) { return; }
2020  if (type==FaceType::Interior)
2021  {
2022  if ( m==L2FaceValues::DoubleValued )
2023  {
2024  DoubleValuedNonconformingTransposeInterpolationInPlace(x);
2025  DoubleValuedConformingAddMultTranspose(x, y);
2026  }
2027  else if ( m==L2FaceValues::SingleValued )
2028  {
2029  SingleValuedNonconformingTransposeInterpolationInPlace(x);
2030  SingleValuedConformingAddMultTranspose(x, y);
2031  }
2032  }
2033  else
2034  {
2035  if ( m==L2FaceValues::DoubleValued )
2036  {
2037  DoubleValuedConformingAddMultTranspose(x, y);
2038  }
2039  else if ( m==L2FaceValues::SingleValued )
2040  {
2041  SingleValuedConformingAddMultTranspose(x, y);
2042  }
2043  }
2044 }
2045 
2046 void NCL2FaceRestriction::FillI(SparseMatrix &mat,
2047  const bool keep_nbr_block) const
2048 {
2049  MFEM_ABORT("Not yet implemented.");
2050 }
2051 
2052 void NCL2FaceRestriction::FillJAndData(const Vector &ea_data,
2053  SparseMatrix &mat,
2054  const bool keep_nbr_block) const
2055 {
2056  MFEM_ABORT("Not yet implemented.");
2057 }
2058 
2059 void NCL2FaceRestriction::AddFaceMatricesToElementMatrices(
2060  const Vector &fea_data,
2061  Vector &ea_data)
2062 const
2063 {
2064  MFEM_ABORT("Not yet implemented.");
2065 }
2066 
2067 int ToLexOrdering(const int dim, const int face_id, const int size1d,
2068  const int index)
2069 {
2070  switch (dim)
2071  {
2072  case 1:
2073  return 0;
2074  case 2:
2075  return ToLexOrdering2D(face_id, size1d, index);
2076  case 3:
2077  return ToLexOrdering3D(face_id, size1d, index%size1d, index/size1d);
2078  default:
2079  MFEM_ABORT("Unsupported dimension.");
2080  return 0;
2081  }
2082 }
2083 
2084 void NCL2FaceRestriction::ComputeScatterIndicesAndOffsets(
2085  const ElementDofOrdering f_ordering,
2086  const FaceType type)
2087 {
2088  Mesh &mesh = *fes.GetMesh();
2089 
2090  // Initialization of the offsets
2091  for (int i = 0; i <= ndofs; ++i)
2092  {
2093  gather_offsets[i] = 0;
2094  }
2095 
2096  // Computation of scatter and offsets indices
2097  int f_ind=0;
2098  for (int f = 0; f < fes.GetNF(); ++f)
2099  {
2100  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
2101  if ( face.IsNonconformingCoarse() )
2102  {
2103  // We skip nonconforming coarse faces as they are treated
2104  // by the corresponding nonconforming fine faces.
2105  continue;
2106  }
2107  else if ( type==FaceType::Interior && face.IsInterior() )
2108  {
2109  SetFaceDofsScatterIndices1(face,f_ind);
2110  if ( m==L2FaceValues::DoubleValued )
2111  {
2112  PermuteAndSetFaceDofsScatterIndices2(face,f_ind);
2113  if ( face.IsConforming() )
2114  {
2115  interpolations.RegisterFaceConformingInterpolation(face,f_ind);
2116  }
2117  else // Non-conforming face
2118  {
2119  interpolations.RegisterFaceCoarseToFineInterpolation(face,f_ind);
2120  }
2121  }
2122  f_ind++;
2123  }
2124  else if ( type==FaceType::Boundary && face.IsBoundary() )
2125  {
2126  SetFaceDofsScatterIndices1(face,f_ind);
2127  if ( m==L2FaceValues::DoubleValued )
2128  {
2129  SetBoundaryDofsScatterIndices2(face,f_ind);
2130  }
2131  f_ind++;
2132  }
2133  }
2134  MFEM_VERIFY(f_ind==nf, "Unexpected number of " <<
2135  (type==FaceType::Interior? "interior" : "boundary") <<
2136  " faces: " << f_ind << " vs " << nf );
2137 
2138  // Summation of the offsets
2139  for (int i = 1; i <= ndofs; ++i)
2140  {
2141  gather_offsets[i] += gather_offsets[i - 1];
2142  }
2143 
2144  // Transform the interpolation matrix map into a contiguous memory structure.
2145  interpolations.LinearizeInterpolatorMapIntoVector();
2146  interpolations.InitializeNCInterpConfig();
2147 }
2148 
2149 void NCL2FaceRestriction::ComputeGatherIndices(
2150  const ElementDofOrdering f_ordering,
2151  const FaceType type)
2152 {
2153  Mesh &mesh = *fes.GetMesh();
2154  // Computation of gather_indices
2155  int f_ind = 0;
2156  for (int f = 0; f < fes.GetNF(); ++f)
2157  {
2158  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
2159  MFEM_ASSERT(!face.IsShared(),
2160  "Unexpected shared face in NCL2FaceRestriction.");
2161  if ( face.IsNonconformingCoarse() )
2162  {
2163  // We skip nonconforming coarse faces as they are treated
2164  // by the corresponding nonconforming fine faces.
2165  continue;
2166  }
2167  else if ( face.IsOfFaceType(type) )
2168  {
2169  SetFaceDofsGatherIndices1(face,f_ind);
2170  if ( m==L2FaceValues::DoubleValued &&
2171  type==FaceType::Interior &&
2172  face.IsInterior() )
2173  {
2174  PermuteAndSetFaceDofsGatherIndices2(face,f_ind);
2175  }
2176  f_ind++;
2177  }
2178  }
2179  MFEM_VERIFY(f_ind==nf, "Unexpected number of " <<
2180  (type==FaceType::Interior? "interior" : "boundary") <<
2181  " faces: " << f_ind << " vs " << nf );
2182 
2183  // Switch back offsets to their correct value
2184  for (int i = ndofs; i > 0; --i)
2185  {
2186  gather_offsets[i] = gather_offsets[i - 1];
2187  }
2188  gather_offsets[0] = 0;
2189 }
2190 
2191 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
void FillI(SparseMatrix &mat) const
Abstract class for all finite elements.
Definition: fe_base.hpp:233
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition: array.hpp:327
Array< NCInterpConfig > nc_interp_config
const Array< NCInterpConfig > & GetNCFaceInterpConfig() const
Return an array containing the interpolation configuration for each face registered with RegisterFace...
int * GetJ()
Definition: table.hpp:114
void forall_2D(int N, int X, int Y, lambda &&body)
Definition: forall.hpp:751
void PermuteAndSetSharedFaceDofsScatterIndices2(const Mesh::FaceInformation &face, const int face_index)
Permute and set the scattering indices of elem2 for the shared face described by the face...
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const override
Add the E-vector degrees of freedom x to the L-vector degrees of freedom y.
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
void SetBoundaryDofsScatterIndices2(const Mesh::FaceInformation &face, const int face_index)
Set the scattering indices of elem2 for the boundary face described by the face.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
std::pair< const DenseMatrix *, int > Key
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
Definition: mesh.cpp:1421
void SetFaceDofsScatterIndices(const Mesh::FaceInformation &face, const int face_index, const ElementDofOrdering f_ordering)
Set the scattering indices of elem1, and increment the offsets for the face described by the face...
MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
Definition: backends.hpp:84
T * GetData()
Returns the data.
Definition: array.hpp:115
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
virtual void DoubleValuedNonconformingMult(const Vector &x, Vector &y) const
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with no...
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:465
Operator that extracts face degrees of freedom for L2 nonconforming spaces.
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
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:3239
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void FillJAndData(const Vector &fea_data, SparseMatrix &mat, const bool keep_nbr_block=false) const
Fill the J and Data arrays of the SparseMatrix corresponding to the sparsity pattern given by this L2...
void SetFaceDofsGatherIndices1(const Mesh::FaceInformation &face, const int face_index)
Set the gathering indices of elem1 for the interior face described by the face.
void AddMultTransposeUnsigned(const Vector &x, Vector &y, const double a=1.0) const override
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector not taking into account signs...
Memory< double > & GetMemoryData()
Definition: sparsemat.hpp:269
const L2FaceValues m
Memory< int > & GetMemoryI()
Definition: sparsemat.hpp:237
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:2841
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const override
Add the E-vector degrees of freedom x to the L-vector degrees of freedom y.
Array< int > gather_indices
void PermuteAndSetFaceDofsGatherIndices2(const Mesh::FaceInformation &face, const int face_index)
Permute and set the gathering indices of elem2 for the interior face described by the face...
const Array< int > & GetDofMap() const
Get an Array<int> that maps lexicographically ordered indices to the indices of the respective nodes/...
Definition: fe_base.hpp:1210
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:3199
L2ElementRestriction(const FiniteElementSpace &)
Memory< int > & GetMemoryJ()
Definition: sparsemat.hpp:253
void DoubleValuedNonconformingInterpolation(Vector &x) const
Apply a change of basis from coarse element basis to fine element basis for the coarse face dofs...
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
Array< InterpConfig > interp_config
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1457
void SingleValuedConformingMult(const Vector &x, Vector &y) const
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with co...
bool IsLocal() const
Return true if the face is a local interior face which is NOT a master nonconforming face...
Definition: mesh.hpp:1683
void SetPointMat(const DenseMatrix &pm)
Set the underlying point matrix describing the transformation.
Definition: eltrans.hpp:401
void SetFaceDofsGatherIndices(const Mesh::FaceInformation &face, const int face_index, const ElementDofOrdering f_ordering)
Set the gathering indices of elem1 for the interior face described by the face.
virtual void AddFaceMatricesToElementMatrices(const Vector &fea_data, Vector &ea_data) const
This methods adds the DG face matrices to the element matrices.
FaceType
Definition: mesh.hpp:45
void RegisterFaceConformingInterpolation(const Mesh::FaceInformation &face, int face_index)
Register the face with face and index face_index as a conforming face for the interpolation of the de...
const FiniteElementSpace & fes
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void PermuteAndSetFaceDofsScatterIndices2(const Mesh::FaceInformation &face, const int face_index)
Permute and set the scattering indices of elem2, and increment the offsets for the face described by ...
Data type sparse matrix.
Definition: sparsemat.hpp:50
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
Operator that extracts Face degrees of freedom for L2 spaces.
ConformingFaceRestriction(const FiniteElementSpace &fes, const ElementDofOrdering f_ordering, const FaceType type, bool build)
Construct a ConformingFaceRestriction.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
void Mult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector.
void LinearizeInterpolatorMapIntoVector()
Transform the interpolation matrix map into a contiguous memory structure.
int * ReadWriteI(bool on_dev=true)
Definition: sparsemat.hpp:243
Bernstein polynomials.
Definition: fe_base.hpp:33
void MultLeftInverse(const Vector &x, Vector &y) const
void MultTransposeUnsigned(const Vector &x, Vector &y) const
Compute MultTranspose without applying signs based on DOF orientations.
virtual void FillI(SparseMatrix &mat, const bool keep_nbr_block=false) const
Fill the I array of SparseMatrix corresponding to the sparsity pattern given by this L2FaceRestrictio...
int * WriteI(bool on_dev=true)
Definition: sparsemat.hpp:241
ElementRestriction(const FiniteElementSpace &, ElementDofOrdering)
Definition: restriction.cpp:27
void Mult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector.
struct mfem::Mesh::FaceInformation::@13 element[2]
const FiniteElementSpace & fes
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition: array.hpp:311
bool IsConforming() const
Return true if the face is a conforming face.
Definition: mesh.hpp:1726
This structure is used as a human readable output format that decipheres the information contained in...
Definition: mesh.hpp:1664
const ElementDofOrdering ordering
int ToLexOrdering(const int dim, const int face_id, const int size1d, const int index)
Convert a dof face index from Native ordering to lexicographic ordering for quads and hexes...
const FiniteElementSpace & fes
Definition: restriction.hpp:45
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
Definition: fespace.hpp:742
Array< int > gather_offsets
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
void MultTranspose(const Vector &x, Vector &y) const override
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition: eltrans.cpp:383
void forall(int N, lambda &&body)
Definition: forall.hpp:742
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
virtual void DoubleValuedConformingMult(const Vector &x, Vector &y) const
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with co...
L2FaceRestriction(const FiniteElementSpace &fes, const ElementDofOrdering f_ordering, const FaceType type, const L2FaceValues m, bool build)
Constructs an L2FaceRestriction.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
void FillJAndData(const Vector &ea_data, SparseMatrix &mat) const
void MultUnsigned(const Vector &x, Vector &y) const
Compute Mult without applying signs based on DOF orientations.
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void SingleValuedConformingAddMultTranspose(const Vector &x, Vector &y) const
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector. Should only be used with con...
void FillSparseMatrix(const Vector &mat_ea, SparseMatrix &mat) const
Fill a Sparse Matrix with Element Matrices.
The ordering method used when the number of unknowns per mesh node (vector dimension) is bigger than ...
Definition: fespace.hpp:29
int * WriteJ(bool on_dev=true)
Definition: sparsemat.hpp:257
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
void CheckFESpace(const ElementDofOrdering f_ordering)
Verify that ConformingFaceRestriction is built from a supported finite element space.
Array< int > scatter_indices1
double a
Definition: lissajous.cpp:41
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:469
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const override
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector.
A standard isoparametric element transformation.
Definition: eltrans.hpp:361
void New(int size)
Allocate host memory for size entries with the current host memory type returned by MemoryManager::Ge...
Operator that extracts face degrees of freedom for H1, ND, or RT FiniteElementSpaces.
Mesh * GetMesh(int type)
Definition: ex29.cpp:218
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
int dim
Definition: ex24.cpp:53
int GetNumInterpolators() const
Return the total number of interpolators.
int FillI(SparseMatrix &mat) const
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
void Mult(const Vector &x, Vector &y) const override
Operator application: y=A(x).
const FiniteElementSpace & fes
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
Lexicographic ordering for tensor-product FiniteElements.
bool IsBoundary() const
Return true if the face is a boundary face.
Definition: mesh.hpp:1706
void DoubleValuedConformingAddMultTranspose(const Vector &x, Vector &y) const
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector. Should only be used with con...
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
RefCoord t[3]
FaceInformation GetFaceInformation(int f) const
Definition: mesh.cpp:1138
ElementConformity conformity
Definition: mesh.hpp:1671
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const override
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector.
Vector data type.
Definition: vector.hpp:58
virtual void GetFaceMap(const int face_id, Array< int > &face_map) const
Return the mapping from lexicographic face DOFs to lexicographic element DOFs for the given local fac...
Definition: fe_base.cpp:494
bool IsShared() const
Return true if the face is a shared interior face which is NOT a master nonconforming face...
Definition: mesh.hpp:1690
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:404
bool IsNonconformingCoarse() const
Return true if the face is a nonconforming coarse face.
Definition: mesh.hpp:1743
int PermuteFaceL2(const int dim, const int face_id1, const int face_id2, const int orientation, const int size1d, const int index)
Compute the dof face index of elem2 corresponding to the given dof face index.
int * HostReadWriteI()
Definition: sparsemat.hpp:249
void SetFaceDofsScatterIndices1(const Mesh::FaceInformation &face, const int face_index)
Set the scattering indices of elem1, and increment the offsets for the face described by the face...
void RegisterFaceCoarseToFineInterpolation(const Mesh::FaceInformation &face, int face_index)
Register the face with face and index face_index as a conforming face for the interpolation of the de...
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
void CheckFESpace(const ElementDofOrdering f_ordering)
Verify that L2FaceRestriction is built from an L2 FESpace.
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:1097
const Vector & GetInterpolators() const
Return an mfem::Vector containing the interpolators in the following format: face_dofs x face_dofs x ...
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
InterpolationManager interpolations
const T & AsConst(const T &a)
Utility function similar to std::as_const in c++17.
Definition: array.hpp:350
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
MemoryType GetMemoryType() const
Return a MemoryType that is currently valid. If both the host and the device pointers are currently v...
Array< int > scatter_indices2
const DenseMatrix * point_matrix
Definition: mesh.hpp:1679
int GetBasisType() const
Definition: fe_base.hpp:1202
VDIM x NQPT x NE (values) / VDIM x DIM x NQPT x NE (grads)
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...
double f(const Vector &p)
double * WriteData(bool on_dev=true)
Definition: sparsemat.hpp:273
NCL2FaceRestriction(const FiniteElementSpace &fes, const ElementDofOrdering f_ordering, const FaceType type, const L2FaceValues m, bool build)
Constructs an NCL2FaceRestriction, this is a specialization of a L2FaceRestriction for nonconforming ...