MFEM  v4.4.0 Finite element discretization library
prestriction.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "../config/config.hpp"
13
14 #ifdef MFEM_USE_MPI
15
16 #include "restriction.hpp"
17 #include "prestriction.hpp"
18 #include "pgridfunc.hpp"
19 #include "pfespace.hpp"
20 #include "fespace.hpp"
21 #include "../general/forall.hpp"
22
23 namespace mfem
24 {
25
27  ElementDofOrdering ordering,
28  FaceType type)
29  : H1FaceRestriction(fes, ordering, type, false),
30  type(type),
31  interpolations(fes, ordering, type)
32 {
33  if (nf==0) { return; }
34  x_interp.UseDevice(true);
35
36  CheckFESpace(ordering);
37
38  ComputeScatterIndicesAndOffsets(ordering, type);
39
40  ComputeGatherIndices(ordering, type);
41 }
42
43 void ParNCH1FaceRestriction::Mult(const Vector &x, Vector &y) const
44 {
45  if (nf==0) { return; }
46  // Assumes all elements have the same number of dofs
47  const int nface_dofs = face_dofs;
48  const int vd = vdim;
49  const bool t = byvdim;
50
51  if ( type==FaceType::Boundary )
52  {
54  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
55  auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
56  MFEM_FORALL(i, nfdofs,
57  {
58  const int dof = i % nface_dofs;
59  const int face = i / nface_dofs;
60  const int idx = d_indices[i];
61  for (int c = 0; c < vd; ++c)
62  {
63  d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
64  }
65  });
66  }
67  else // type==FaceType::Interior
68  {
70  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
71  auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
74  const int nc_size = interpolations.GetNumInterpolators();
75  auto d_interp = Reshape(interpolators, nface_dofs, nface_dofs, nc_size);
76  static constexpr int max_nd = 1024;
77  MFEM_VERIFY(nface_dofs<=max_nd, "Too many degrees of freedom.");
78  MFEM_FORALL_3D(face, nf, nface_dofs, 1, 1,
79  {
80  MFEM_SHARED double dof_values[max_nd];
81  const InterpConfig conf = interp_config_ptr[face];
82  const int master_side = conf.master_side;
83  const int interp_index = conf.index;
84  const int side = 0;
85  if ( !conf.is_non_conforming || side!=master_side )
86  {
88  {
89  const int i = face*nface_dofs + dof;
90  const int idx = d_indices[i];
91  for (int c = 0; c < vd; ++c)
92  {
93  d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
94  }
95  }
96  }
97  else // Interpolation from coarse to fine
98  {
99  for (int c = 0; c < vd; ++c)
100  {
101  // Load the face dofs in shared memory
103  {
104  const int i = face*nface_dofs + dof;
105  const int idx = d_indices[i];
106  dof_values[dof] = d_x(t?c:idx, t?idx:c);
107  }
109  // Apply the interpolation to the face dofs
111  {
112  double res = 0.0;
113  for (int dof_in = 0; dof_in<nface_dofs; dof_in++)
114  {
115  res += d_interp(dof_out, dof_in, interp_index)*
116  dof_values[dof_in];
117  }
118  d_y(dof_out, c, face) = res;
119  }
121  }
122  }
123  });
124  }
125 }
126
127 void ParNCH1FaceRestriction::AddMultTranspose(const Vector &x, Vector &y) const
128 {
129  if (nf==0) { return; }
130  if (x_interp.Size()==0)
131  {
132  x_interp.SetSize(x.Size());
133  }
134  x_interp = x;
135  // Assumes all elements have the same number of dofs
136  const int nface_dofs = face_dofs;
137  const int vd = vdim;
138  const bool t = byvdim;
139  if ( type==FaceType::Interior )
140  {
141  // Interpolation from slave to master face dofs
142  auto d_x = Reshape(x_interp.ReadWrite(), nface_dofs, vd, nf);
145  const int nc_size = interpolations.GetNumInterpolators();
146  auto d_interp = Reshape(interpolators, nface_dofs, nface_dofs, nc_size);
147  static constexpr int max_nd = 1024;
148  MFEM_VERIFY(nface_dofs<=max_nd, "Too many degrees of freedom.");
149  MFEM_FORALL_3D(face, nf, nface_dofs, 1, 1,
150  {
151  MFEM_SHARED double dof_values[max_nd];
152  const InterpConfig conf = interp_config_ptr[face];
153  const int master_side = conf.master_side;
154  const int interp_index = conf.index;
155  if ( conf.is_non_conforming && master_side==0 )
156  {
157  // Interpolation from fine to coarse
158  for (int c = 0; c < vd; ++c)
159  {
161  {
162  dof_values[dof] = d_x(dof, c, face);
163  }
166  {
167  double res = 0.0;
168  for (int dof_in = 0; dof_in<nface_dofs; dof_in++)
169  {
170  res += d_interp(dof_in, dof_out, interp_index)*dof_values[dof_in];
171  }
172  d_x(dof_out, c, face) = res;
173  }
175  }
176  }
177  });
178  }
179
180  // Gathering of face dofs into element dofs
183  auto d_x = Reshape(x_interp.Read(), nface_dofs, vd, nf);
184  auto d_y = Reshape(y.ReadWrite(), t?vd:ndofs, t?ndofs:vd);
185  MFEM_FORALL(i, ndofs,
186  {
187  const int offset = d_offsets[i];
188  const int next_offset = d_offsets[i + 1];
189  for (int c = 0; c < vd; ++c)
190  {
191  double dof_value = 0;
192  for (int j = offset; j < next_offset; ++j)
193  {
194  int idx_j = d_indices[j];
195  dof_value += d_x(idx_j % nface_dofs, c, idx_j / nface_dofs);
196  }
197  d_y(t?c:i,t?i:c) += dof_value;
198  }
199  });
200 }
201
202 void ParNCH1FaceRestriction::ComputeScatterIndicesAndOffsets(
203  const ElementDofOrdering ordering,
204  const FaceType face_type)
205 {
206  Mesh &mesh = *fes.GetMesh();
207
208  // Initialization of the offsets
209  for (int i = 0; i <= ndofs; ++i)
210  {
211  gather_offsets[i] = 0;
212  }
213
214  // Computation of scatter indices and offsets
215  int f_ind = 0;
216  for (int f = 0; f < mesh.GetNumFacesWithGhost(); ++f)
217  {
218  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
219  if ( face.IsNonconformingCoarse() )
220  {
221  // We skip nonconforming coarse faces as they are treated
222  // by the corresponding nonconforming fine faces.
223  continue;
224  }
225  else if (face_type==FaceType::Interior && face.IsInterior())
226  {
227  if ( face.IsConforming() )
228  {
229  interpolations.RegisterFaceConformingInterpolation(face,f_ind);
230  SetFaceDofsScatterIndices(face, f_ind, ordering);
231  f_ind++;
232  }
233  else // Non-conforming face
234  {
235  SetFaceDofsScatterIndices(face, f_ind, ordering);
236  if ( face.element[0].conformity==Mesh::ElementConformity::Superset )
237  {
238  // In this case the local face is the master (coarse) face, thus
239  // we need to interpolate the values on the slave (fine) face.
240  interpolations.RegisterFaceCoarseToFineInterpolation(face,f_ind);
241  }
242  else
243  {
244  // Treated as a conforming face since we only extract values from
245  // the local slave (fine) face.
246  interpolations.RegisterFaceConformingInterpolation(face,f_ind);
247  }
248  f_ind++;
249  }
250  }
251  else if (face_type==FaceType::Boundary && face.IsBoundary())
252  {
253  SetFaceDofsScatterIndices(face, f_ind, ordering);
254  f_ind++;
255  }
256  }
257  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
258
259  // Summation of the offsets
260  for (int i = 1; i <= ndofs; ++i)
261  {
262  gather_offsets[i] += gather_offsets[i - 1];
263  }
264
265  // Transform the interpolation matrix map into a contiguous memory structure.
266  interpolations.LinearizeInterpolatorMapIntoVector();
267 }
268
269 void ParNCH1FaceRestriction::ComputeGatherIndices(
270  const ElementDofOrdering ordering,
271  const FaceType face_type)
272 {
273  Mesh &mesh = *fes.GetMesh();
274
275  // Computation of gather_indices
276  int f_ind = 0;
277  for (int f = 0; f < mesh.GetNumFacesWithGhost(); ++f)
278  {
279  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
280  if ( face.IsNonconformingCoarse() )
281  {
282  // We skip nonconforming coarse faces as they are treated
283  // by the corresponding nonconforming fine faces.
284  continue;
285  }
286  else if (face.IsOfFaceType(face_type))
287  {
288  SetFaceDofsGatherIndices(face, f_ind, ordering);
289  f_ind++;
290  }
291  }
292  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
293
294  // Reset offsets to their correct value
295  for (int i = ndofs; i > 0; --i)
296  {
297  gather_offsets[i] = gather_offsets[i - 1];
298  }
299  gather_offsets[0] = 0;
300 }
301
302 ParL2FaceRestriction::ParL2FaceRestriction(const ParFiniteElementSpace &fes,
303  ElementDofOrdering ordering,
304  FaceType type,
305  L2FaceValues m,
306  bool build)
307  : L2FaceRestriction(fes, ordering, type, m, false)
308 {
309  if (!build) { return; }
310  if (nf==0) { return; }
311
312  CheckFESpace(ordering);
313
314  ComputeScatterIndicesAndOffsets(ordering, type);
315
316  ComputeGatherIndices(ordering, type);
317 }
318
320  ElementDofOrdering ordering,
321  FaceType type,
322  L2FaceValues m)
323  : ParL2FaceRestriction(fes, ordering, type, m, true)
324 { }
325
327  const Vector& x, Vector& y) const
328 {
329  MFEM_ASSERT(
331  "This method should be called when m == L2FaceValues::DoubleValued.");
332  const ParFiniteElementSpace &pfes =
333  static_cast<const ParFiniteElementSpace&>(this->fes);
334  ParGridFunction x_gf;
335  x_gf.MakeRef(const_cast<ParFiniteElementSpace*>(&pfes),
336  const_cast<Vector&>(x), 0);
337  x_gf.ExchangeFaceNbrData();
338
339  // Assumes all elements have the same number of dofs
340  const int nface_dofs = face_dofs;
341  const int vd = vdim;
342  const bool t = byvdim;
343  const int threshold = ndofs;
344  const int nsdofs = pfes.GetFaceNbrVSize();
347  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
349  t?vd:nsdofs, t?nsdofs:vd);
350  auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
351  MFEM_FORALL(i, nfdofs,
352  {
353  const int dof = i % nface_dofs;
354  const int face = i / nface_dofs;
355  const int idx1 = d_indices1[i];
356  for (int c = 0; c < vd; ++c)
357  {
358  d_y(dof, c, 0, face) = d_x(t?c:idx1, t?idx1:c);
359  }
360  const int idx2 = d_indices2[i];
361  for (int c = 0; c < vd; ++c)
362  {
363  if (idx2>-1 && idx2<threshold) // interior face
364  {
365  d_y(dof, c, 1, face) = d_x(t?c:idx2, t?idx2:c);
366  }
367  else if (idx2>=threshold) // shared boundary
368  {
369  d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
370  t?(idx2-threshold):c);
371  }
372  else // true boundary
373  {
374  d_y(dof, c, 1, face) = 0.0;
375  }
376  }
377  });
378 }
379
380 void ParL2FaceRestriction::Mult(const Vector& x, Vector& y) const
381 {
382  if (nf==0) { return; }
384  {
386  }
387  else
388  {
390  }
391 }
392
393 static MFEM_HOST_DEVICE int AddNnz(const int iE, int *I, const int dofs)
394 {
396  return val;
397 }
398
400  const bool keep_nbr_block) const
401 {
402  if (keep_nbr_block)
403  {
404  return L2FaceRestriction::FillI(mat, keep_nbr_block);
405  }
406  const int nface_dofs = face_dofs;
407  const int Ndofs = ndofs;
411  MFEM_FORALL(fdof, nf*nface_dofs,
412  {
413  const int f = fdof/nface_dofs;
414  const int iF = fdof%nface_dofs;
415  const int iE1 = d_indices1[f*nface_dofs+iF];
416  if (iE1 < Ndofs)
417  {
419  }
420  const int iE2 = d_indices2[f*nface_dofs+iF];
421  if (iE2 < Ndofs)
422  {
424  }
425  });
426 }
427
429  SparseMatrix &face_mat) const
430 {
431  const int nface_dofs = face_dofs;
432  const int Ndofs = ndofs;
437  MFEM_FORALL(i, ne*elem_dofs*vdim+1,
438  {
439  I_face[i] = 0;
440  });
441  MFEM_FORALL(fdof, nf*nface_dofs,
442  {
443  const int f = fdof/nface_dofs;
444  const int iF = fdof%nface_dofs;
445  const int iE1 = d_indices1[f*nface_dofs+iF];
446  if (iE1 < Ndofs)
447  {
448  for (int jF = 0; jF < nface_dofs; jF++)
449  {
450  const int jE2 = d_indices2[f*nface_dofs+jF];
451  if (jE2 < Ndofs)
452  {
454  }
455  else
456  {
458  }
459  }
460  }
461  const int iE2 = d_indices2[f*nface_dofs+iF];
462  if (iE2 < Ndofs)
463  {
464  for (int jF = 0; jF < nface_dofs; jF++)
465  {
466  const int jE1 = d_indices1[f*nface_dofs+jF];
467  if (jE1 < Ndofs)
468  {
470  }
471  else
472  {
474  }
475  }
476  }
477  });
478 }
479
481  SparseMatrix &mat,
482  const bool keep_nbr_block) const
483 {
484  if (keep_nbr_block)
485  {
486  return L2FaceRestriction::FillJAndData(ea_data, mat, keep_nbr_block);
487  }
488  const int nface_dofs = face_dofs;
489  const int Ndofs = ndofs;
492  auto mat_fea = Reshape(ea_data.Read(), nface_dofs, nface_dofs, 2, nf);
494  auto J = mat.WriteJ();
495  auto Data = mat.WriteData();
496  MFEM_FORALL(fdof, nf*nface_dofs,
497  {
498  const int f = fdof/nface_dofs;
499  const int iF = fdof%nface_dofs;
500  const int iE1 = d_indices1[f*nface_dofs+iF];
501  if (iE1 < Ndofs)
502  {
503  const int offset = AddNnz(iE1,I,nface_dofs);
504  for (int jF = 0; jF < nface_dofs; jF++)
505  {
506  const int jE2 = d_indices2[f*nface_dofs+jF];
507  J[offset+jF] = jE2;
508  Data[offset+jF] = mat_fea(jF,iF,1,f);
509  }
510  }
511  const int iE2 = d_indices2[f*nface_dofs+iF];
512  if (iE2 < Ndofs)
513  {
514  const int offset = AddNnz(iE2,I,nface_dofs);
515  for (int jF = 0; jF < nface_dofs; jF++)
516  {
517  const int jE1 = d_indices1[f*nface_dofs+jF];
518  J[offset+jF] = jE1;
519  Data[offset+jF] = mat_fea(jF,iF,0,f);
520  }
521  }
522  });
523 }
524
526  SparseMatrix &mat,
527  SparseMatrix &face_mat) const
528 {
529  const int nface_dofs = face_dofs;
530  const int Ndofs = ndofs;
533  auto mat_fea = Reshape(ea_data.Read(), nface_dofs, nface_dofs, 2, nf);
536  auto J = mat.WriteJ();
537  auto J_face = face_mat.WriteJ();
538  auto Data = mat.WriteData();
539  auto Data_face = face_mat.WriteData();
540  MFEM_FORALL(fdof, nf*nface_dofs,
541  {
542  const int f = fdof/nface_dofs;
543  const int iF = fdof%nface_dofs;
544  const int iE1 = d_indices1[f*nface_dofs+iF];
545  if (iE1 < Ndofs)
546  {
547  for (int jF = 0; jF < nface_dofs; jF++)
548  {
549  const int jE2 = d_indices2[f*nface_dofs+jF];
550  if (jE2 < Ndofs)
551  {
552  const int offset = AddNnz(iE1,I,1);
553  J[offset] = jE2;
554  Data[offset] = mat_fea(jF,iF,1,f);
555  }
556  else
557  {
558  const int offset = AddNnz(iE1,I_face,1);
559  J_face[offset] = jE2-Ndofs;
560  Data_face[offset] = mat_fea(jF,iF,1,f);
561  }
562  }
563  }
564  const int iE2 = d_indices2[f*nface_dofs+iF];
565  if (iE2 < Ndofs)
566  {
567  for (int jF = 0; jF < nface_dofs; jF++)
568  {
569  const int jE1 = d_indices1[f*nface_dofs+jF];
570  if (jE1 < Ndofs)
571  {
572  const int offset = AddNnz(iE2,I,1);
573  J[offset] = jE1;
574  Data[offset] = mat_fea(jF,iF,0,f);
575  }
576  else
577  {
578  const int offset = AddNnz(iE2,I_face,1);
579  J_face[offset] = jE1-Ndofs;
580  Data_face[offset] = mat_fea(jF,iF,0,f);
581  }
582  }
583  }
584  });
585 }
586
587 void ParL2FaceRestriction::ComputeScatterIndicesAndOffsets(
588  const ElementDofOrdering ordering,
589  const FaceType type)
590 {
591  Mesh &mesh = *fes.GetMesh();
592  const ParFiniteElementSpace &pfes =
593  static_cast<const ParFiniteElementSpace&>(this->fes);
594
595  // Initialization of the offsets
596  for (int i = 0; i <= ndofs; ++i)
597  {
598  gather_offsets[i] = 0;
599  }
600
601  // Computation of scatter indices and offsets
602  int f_ind=0;
603  for (int f = 0; f < pfes.GetNF(); ++f)
604  {
605  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
606  if (type==FaceType::Interior && face.IsInterior())
607  {
608  SetFaceDofsScatterIndices1(face,f_ind);
610  {
611  if (face.IsShared())
612  {
614  }
615  else
616  {
618  }
619  }
620  f_ind++;
621  }
622  else if (type==FaceType::Boundary && face.IsBoundary())
623  {
624  SetFaceDofsScatterIndices1(face,f_ind);
626  {
627  SetBoundaryDofsScatterIndices2(face,f_ind);
628  }
629  f_ind++;
630  }
631  }
632  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
633
634  // Summation of the offsets
635  for (int i = 1; i <= ndofs; ++i)
636  {
637  gather_offsets[i] += gather_offsets[i - 1];
638  }
639 }
640
641
642 void ParL2FaceRestriction::ComputeGatherIndices(
643  const ElementDofOrdering ordering,
644  const FaceType type)
645 {
646  Mesh &mesh = *fes.GetMesh();
647
648  // Computation of gather_indices
649  int f_ind = 0;
650  for (int f = 0; f < fes.GetNF(); ++f)
651  {
652  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
653  if (face.IsOfFaceType(type))
654  {
655  SetFaceDofsGatherIndices1(face,f_ind);
657  type==FaceType::Interior &&
658  face.IsLocal())
659  {
661  }
662  f_ind++;
663  }
664  }
665  MFEM_VERIFY(f_ind==nf, "Unexpected number of faces.");
666
667  // Reset offsets to their correct value
668  for (int i = ndofs; i > 0; --i)
669  {
670  gather_offsets[i] = gather_offsets[i - 1];
671  }
672  gather_offsets[0] = 0;
673 }
674
676  ElementDofOrdering ordering,
677  FaceType type,
678  L2FaceValues m)
679  : L2FaceRestriction(fes, ordering, type, m, false),
680  NCL2FaceRestriction(fes, ordering, type, m, false),
681  ParL2FaceRestriction(fes, ordering, type, m, false)
682 {
683  if (nf==0) { return; }
684  x_interp.UseDevice(true);
685
686  CheckFESpace(ordering);
687
688  ComputeScatterIndicesAndOffsets(ordering, type);
689
690  ComputeGatherIndices(ordering, type);
691 }
692
694  const Vector& x, Vector& y) const
695 {
696  MFEM_ASSERT(
698  "This method should be called when m == L2FaceValues::SingleValued.");
699  // Assumes all elements have the same number of dofs
700  const int nface_dofs = face_dofs;
701  const int vd = vdim;
702  const bool t = byvdim;
703  const int threshold = ndofs;
705  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
706  auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
709  const int nc_size = interpolations.GetNumInterpolators();
710  auto d_interp = Reshape(interpolators, nface_dofs, nface_dofs, nc_size);
711  static constexpr int max_nd = 16*16;
712  MFEM_VERIFY(nface_dofs<=max_nd, "Too many degrees of freedom.");
713  MFEM_FORALL_3D(face, nf, nface_dofs, 1, 1,
714  {
715  MFEM_SHARED double dof_values[max_nd];
716  const InterpConfig conf = interp_config_ptr[face];
717  const int master_side = conf.master_side;
718  const int interp_index = conf.index;
719  const int side = 0;
720  if ( !conf.is_non_conforming || side!=master_side )
721  {
723  {
724  const int i = face*nface_dofs + dof;
725  const int idx = d_indices1[i];
726  if (idx>-1 && idx<threshold) // interior face
727  {
728  for (int c = 0; c < vd; ++c)
729  {
730  d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
731  }
732  }
733  else // true boundary
734  {
735  for (int c = 0; c < vd; ++c)
736  {
737  d_y(dof, c, face) = 0.0;
738  }
739  }
740  }
741  }
742  else // Interpolation from coarse to fine
743  {
744  for (int c = 0; c < vd; ++c)
745  {
747  {
748  const int i = face*nface_dofs + dof;
749  const int idx = d_indices1[i];
750  if (idx>-1 && idx<threshold) // interior face
751  {
752  dof_values[dof] = d_x(t?c:idx, t?idx:c);
753  }
754  else // true boundary
755  {
756  dof_values[dof] = 0.0;
757  }
758  }
761  {
762  double res = 0.0;
763  for (int dof_in = 0; dof_in<nface_dofs; dof_in++)
764  {
765  res += d_interp(dof_out, dof_in, interp_index)*dof_values[dof_in];
766  }
767  d_y(dof_out, c, face) = res;
768  }
770  }
771  }
772  });
773 }
774
776  const Vector& x, Vector& y) const
777 {
778  MFEM_ASSERT(
780  "This method should be called when m == L2FaceValues::DoubleValued.");
781  const ParFiniteElementSpace &pfes =
782  static_cast<const ParFiniteElementSpace&>(this->fes);
783  ParGridFunction x_gf;
784  x_gf.MakeRef(const_cast<ParFiniteElementSpace*>(&pfes),
785  const_cast<Vector&>(x), 0);
786  x_gf.ExchangeFaceNbrData();
787
788  // Assumes all elements have the same number of dofs
789  const int nface_dofs = face_dofs;
790  const int vd = vdim;
791  const bool t = byvdim;
792  const int threshold = ndofs;
793  const int nsdofs = pfes.GetFaceNbrVSize();
796  auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
798  t?vd:nsdofs, t?nsdofs:vd);
799  auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
802  const int nc_size = interpolations.GetNumInterpolators();
803  auto d_interp = Reshape(interpolators, nface_dofs, nface_dofs, nc_size);
804  static constexpr int max_nd = 1024;
805  MFEM_VERIFY(nface_dofs<=max_nd, "Too many degrees of freedom.");
806  MFEM_FORALL_3D(face, nf, nface_dofs, 1, 1,
807  {
808  MFEM_SHARED double dof_values[max_nd];
809  const InterpConfig conf = interp_config_ptr[face];
810  const int master_side = conf.master_side;
811  const int interp_index = conf.index;
812  for (int side = 0; side < 2; side++)
813  {
814  if ( !conf.is_non_conforming || side!=master_side )
815  {
816  // No interpolation
818  {
819  const int i = face*nface_dofs + dof;
820  const int idx = side==0 ? d_indices1[i] : d_indices2[i];
821  if (idx>-1 && idx<threshold) // local interior face
822  {
823  for (int c = 0; c < vd; ++c)
824  {
825  d_y(dof, c, side, face) = d_x(t?c:idx, t?idx:c);
826  }
827  }
828  else if (idx>=threshold) // shared interior face
829  {
830  const int sidx = idx-threshold;
831  for (int c = 0; c < vd; ++c)
832  {
833  d_y(dof, c, side, face) = d_x_shared(t?c:sidx, t?sidx:c);
834  }
835  }
836  else // true boundary
837  {
838  for (int c = 0; c < vd; ++c)
839  {
840  d_y(dof, c, side, face) = 0.0;
841  }
842  }
843  }
844  }
845  else // Interpolation from coarse to fine
846  {
847  for (int c = 0; c < vd; ++c)
848  {
850  {
851  const int i = face*nface_dofs + dof;
852  const int idx = side==0 ? d_indices1[i] : d_indices2[i];
853  if (idx>-1 && idx<threshold) // local interior face
854  {
855  dof_values[dof] = d_x(t?c:idx, t?idx:c);
856  }
857  else if (idx>=threshold) // shared interior face
858  {
859  const int sidx = idx-threshold;
860  dof_values[dof] = d_x_shared(t?c:sidx, t?sidx:c);
861  }
862  else // true boundary
863  {
864  dof_values[dof] = 0.0;
865  }
866  }
869  {
870  double res = 0.0;
871  for (int dof_in = 0; dof_in<nface_dofs; dof_in++)
872  {
873  res += d_interp(dof_out, dof_in, interp_index)*dof_values[dof_in];
874  }
875  d_y(dof_out, c, side, face) = res;
876  }
878  }
879  }
880  }
881  });
882 }
883
885 {
886  if (nf==0) { return; }
888  {
890  }
891  else if ( type==FaceType::Boundary && m==L2FaceValues::DoubleValued )
892  {
894  }
895  else if ( type==FaceType::Interior && m==L2FaceValues::SingleValued )
896  {
898  }
899  else if ( type==FaceType::Boundary && m==L2FaceValues::SingleValued )
900  {
902  }
903  else
904  {
905  MFEM_ABORT("Unknown type and multiplicity combination.");
906  }
907 }
908
910 {
911  if (nf==0) { return; }
912  if (type==FaceType::Interior)
913  {
915  {
918  }
919  else // Single Valued
920  {
923  }
924  }
925  else
926  {
928  {
930  }
931  else // Single valued
932  {
934  }
935  }
936 }
937
939  const bool keep_nbr_block) const
940 {
941  MFEM_ABORT("Not yet implemented.");
942 }
943
945  SparseMatrix &face_mat) const
946 {
947  MFEM_ABORT("Not yet implemented.");
948 }
949
951  SparseMatrix &mat,
952  const bool keep_nbr_block) const
953 {
954  MFEM_ABORT("Not yet implemented.");
955 }
956
958  SparseMatrix &mat,
959  SparseMatrix &face_mat) const
960 {
961  MFEM_ABORT("Not yet implemented.");
962 }
963
964 void ParNCL2FaceRestriction::ComputeScatterIndicesAndOffsets(
965  const ElementDofOrdering ordering,
966  const FaceType type)
967 {
968  Mesh &mesh = *fes.GetMesh();
969
970  // Initialization of the offsets
971  for (int i = 0; i <= ndofs; ++i)
972  {
973  gather_offsets[i] = 0;
974  }
975
976  // Computation of scatter and offsets indices
977  int f_ind=0;
978  for (int f = 0; f < mesh.GetNumFacesWithGhost(); ++f)
979  {
980  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
981  if ( face.IsNonconformingCoarse() )
982  {
983  // We skip nonconforming coarse faces as they are treated
984  // by the corresponding nonconforming fine faces.
985  continue;
986  }
987  else if ( type==FaceType::Interior && face.IsInterior() )
988  {
989  if ( face.IsConforming() )
990  {
992  SetFaceDofsScatterIndices1(face,f_ind);
994  {
995  if ( face.IsShared() )
996  {
998  }
999  else
1000  {
1002  }
1003  }
1004  }
1005  else // Non-conforming face
1006  {
1008  SetFaceDofsScatterIndices1(face,f_ind);
1010  {
1011  if ( face.IsShared() )
1012  {
1014  }
1015  else // local nonconforming slave
1016  {
1018  }
1019  }
1020  }
1021  f_ind++;
1022  }
1023  else if (type==FaceType::Boundary && face.IsBoundary())
1024  {
1025  SetFaceDofsScatterIndices1(face,f_ind);
1027  {
1028  SetBoundaryDofsScatterIndices2(face,f_ind);
1029  }
1030  f_ind++;
1031  }
1032  }
1033  MFEM_VERIFY(f_ind==nf, "Unexpected number of " <<
1034  (type==FaceType::Interior? "interior" : "boundary") <<
1035  " faces: " << f_ind << " vs " << nf );
1036
1037  // Summation of the offsets
1038  for (int i = 1; i <= ndofs; ++i)
1039  {
1040  gather_offsets[i] += gather_offsets[i - 1];
1041  }
1042
1043  // Transform the interpolation matrix map into a contiguous memory structure.
1045 }
1046
1047 void ParNCL2FaceRestriction::ComputeGatherIndices(
1048  const ElementDofOrdering ordering,
1049  const FaceType type)
1050 {
1051  Mesh &mesh = *fes.GetMesh();
1052
1053  // Computation of gather_indices
1054  int f_ind = 0;
1055  for (int f = 0; f < mesh.GetNumFacesWithGhost(); ++f)
1056  {
1057  Mesh::FaceInformation face = mesh.GetFaceInformation(f);
1058  if ( face.IsNonconformingCoarse() )
1059  {
1060  // We skip nonconforming coarse faces as they are treated
1061  // by the corresponding nonconforming fine faces.
1062  continue;
1063  }
1064  else if ( face.IsOfFaceType(type) )
1065  {
1066  SetFaceDofsGatherIndices1(face,f_ind);
1068  type==FaceType::Interior &&
1069  face.IsLocal())
1070  {
1072  }
1073  f_ind++;
1074  }
1075  }
1076  MFEM_VERIFY(f_ind==nf, "Unexpected number of " <<
1077  (type==FaceType::Interior? "interior" : "boundary") <<
1078  " faces: " << f_ind << " vs " << nf );
1079
1080  // Switch back offsets to their correct value
1081  for (int i = ndofs; i > 0; --i)
1082  {
1083  gather_offsets[i] = gather_offsets[i - 1];
1084  }
1085  gather_offsets[0] = 0;
1086 }
1087
1088 } // namespace mfem
1089
1090 #endif
void Mult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector.
FaceInformation GetFaceInformation(int f) const
Definition: mesh.cpp:1130
Operator that extracts Face degrees of freedom in parallel.
void DoubleValuedNonconformingMult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with no...
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...
Operator that extracts Face degrees of freedom for H1 FiniteElementSpaces.
void Mult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector.
ParL2FaceRestriction(const ParFiniteElementSpace &fes, ElementDofOrdering ordering, FaceType type, L2FaceValues m, bool build)
Constructs an ParL2FaceRestriction.
ParNCH1FaceRestriction(const ParFiniteElementSpace &fes, ElementDofOrdering ordering, FaceType type)
Constructs an ParNCH1FaceRestriction.
void SetBoundaryDofsScatterIndices2(const Mesh::FaceInformation &face, const int face_index)
Set the scattering indices of elem2 for the boundary face described by the face.
Definition: backends.hpp:84
InterpolationManager interpolations
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Operator that extracts face degrees of freedom for L2 nonconforming spaces.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
void SetFaceDofsGatherIndices1(const Mesh::FaceInformation &face, const int face_index)
Set the gathering indices of elem1 for the interior face described by the face.
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:117
const L2FaceValues m
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 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...
ParNCL2FaceRestriction(const ParFiniteElementSpace &fes, ElementDofOrdering ordering, FaceType type, L2FaceValues m=L2FaceValues::DoubleValued)
Constructs an ParNCL2FaceRestriction.
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...
DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims...dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:102
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
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
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:46
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...
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:433
Operator that extracts Face degrees of freedom for L2 spaces.
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:446
void CheckFESpace(const ElementDofOrdering ordering)
Verify that H1FaceRestriction is build from an H1 FESpace.
void LinearizeInterpolatorMapIntoVector()
Transform the interpolation matrix map into a contiguous memory structure.
Definition: sparsemat.hpp:230
const T * Read(bool on_dev=true) const
Definition: array.hpp:304
void FillJAndData(const Vector &fea_data, SparseMatrix &mat, SparseMatrix &face_mat) const
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
Definition: fespace.hpp:596
Vector & FaceNbrData()
Definition: pgridfunc.hpp:208
Array< int > gather_offsets
void DoubleValuedConformingMult(const Vector &x, Vector &y) const override
Scatter the degrees of freedom, i.e. goes from L-Vector to face E-Vector. Should only be used with co...
void SingleValuedNonconformingMult(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...
int GetNumInterpolators() const
Return the total number of interpolators.
int * WriteJ(bool on_dev=true)
Definition: sparsemat.hpp:244
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...
Array< int > scatter_indices1
void CheckFESpace(const ElementDofOrdering ordering)
Verify that L2FaceRestriction is build from an L2 FESpace.
int GetNumFacesWithGhost() const
Return the number of faces (3D), edges (2D) or vertices (1D) including ghost faces.
Definition: mesh.cpp:5356
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:66
int GetFaceNbrVSize() const
Definition: pfespace.hpp:394
RefCoord t[3]
void SingleValuedNonconformingTransposeInterpolation(const Vector &x) const
Apply a change of basis from fine element basis to coarse element basis for the coarse face dofs...
Vector data type.
Definition: vector.hpp:60
void DoubleValuedNonconformingTransposeInterpolation(const Vector &x) const
Apply a change of basis from fine element basis to coarse element basis for the coarse face dofs...
void FillJAndData(const Vector &ea_data, SparseMatrix &mat, SparseMatrix &face_mat) const
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 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...
Class for parallel grid function.
Definition: pgridfunc.hpp:32
uint32_t is_non_conforming
const Array< InterpConfig > & GetFaceInterpConfig() const
Return an array containing the interpolation configuration for each face registered with RegisterFace...
void AddMultTranspose(const Vector &x, Vector &y) const override
Gather the degrees of freedom, i.e. goes from face E-Vector to L-Vector.
virtual const double * Read(bool on_dev=true) const
Definition: vector.hpp:438
Array< int > scatter_indices
InterpolationManager interpolations
void FillI(SparseMatrix &mat, const bool keep_nbr_block=false) const override
Array< int > scatter_indices2
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...
double f(const Vector &p)
double * WriteData(bool on_dev=true)
Definition: sparsemat.hpp:260
void FillI(SparseMatrix &mat, const bool keep_nbr_block=false) const override
Fill the I array of SparseMatrix corresponding to the sparsity pattern given by this ParNCL2FaceRestr...
const Vector & GetInterpolators() const
Return an mfem::Vector containing the interpolators in the following format: face_dofs x face_dofs x ...